Science, Faculty of
DSpace
UBCV
Salehi, Sohrab
Doctor of Philosophy - PhD
University of British Columbia
Tumour fitness landscapes underpin selection in cancer, impacting evolution and response to treatment. Quantitative fitness modelling of cancer cells has numerous and diverse implications: attributing clonal dynamics to drift or selection, identifying the determinants of clonal expansion, and forecasting tumour growth trajectories. Why and how drug resistance evolves is among the key unresolved areas of investigation that require advanced understanding of fitness in cancer. Longitudinal xenoengraftment interrogated via next generation single cell sequencing (SCS) has enabled more accurate, quantitative measurements of tumours as they evolve. This process generates timestamped samples comprising thousands of cells each measured at thousands of copy number aberrations (CNA). Major analytical challenges introduced by this new datatype include (i) how to identify biologically meaningful groups of cells (i.e., clones) across multiple timepoints, and (ii) how to quantitatively reason about the underlying evolutionary forces acting on the clones via their observed dynamics. To address the first problem, we describe and provide supplementary tools for sitka, a scalable Bayesian phylogenetic inference method in Chapter 2. It resolves the clonal structure of a heterogeneous tumour cell population sampled over multiple timepoints by reconstructing the evolutionary relationship between single cells from their inferred CNA profiles. We then develop Lumberjack, a tree-cutting algorithm, and use it to assign cells to clones. We address the second problem in Chapter 3 by developing fitClone, a Bayesian probabilistic framework that ascribes quantitative selection coefficients to individual cancer clones and forecasts competitive clonal dynamics over time. In Chapter 4, we exemplify the computational models introduced above on real-world data collected from cancer cells over a multi-year period to verify two key hypotheses, that (i) clonal dynamics in a pre-treatment triple negative breast cancer (TNBC) tumour is quantifiably reproducible, and that (ii) the fitness landscape is reversed under early response to cisplatin treatment. Our results show that population genetic modelling of timeseries tumour measurements to predict clonal evolution is tractable. Further study with timeseries modelling will provide insight into therapeutic strategies promoting early intervention, drug combinations and evolution-aware approaches to clinical management.
Quantitative fitness modelling incancer using single cell timeseriespopulation dynamicsbySohrab SalehiB.Sc., Amirkabir University of Technology, 2011M.Sc., The University of British Columbia, 2015A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Bioinformatics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)February 2021© Sohrab Salehi 2021The following individuals certify that they have read, and recommend to the Faculty of Graduateand Postdoctoral Studies for acceptance, the thesis entitled:Quantitative fitness modelling in cancer using single cell timeseries population dynam-icssubmitted by Sohrab Salehi in partial fulfillment of the requirements for the degree of Doctorof Philosophy in Bioinformatics.Examining Committee:Sohrab Shah, Department of Pathology and Laboratory Medicine, UBCCo-supervisorAlexandre Bouchard-Coˆte´, Department of Statistics, UBCCo-supervisorSamuel Aparicio, Department of Molecular Oncology, UBCSupervisory Committee MemberRyan Brinkman, Department of Medical Genetics, UBCUniversity ExaminerGeoffrey Schiebinger, Department of Mathematics, UBCUniversity ExaminerAdditional Supervisory Committee Members:Sarah Otto, Department of Zoology, UBCSupervisory Committee MemberiiAbstractTumour fitness landscapes underpin selection in cancer, impacting evolution and response to treat-ment. Quantitative fitness modelling of cancer cells has numerous and diverse implications: at-tributing clonal dynamics to drift or selection, identifying the determinants of clonal expansion,and forecasting tumour growth trajectories. Why and how drug resistance evolves is among thekey unresolved areas of investigation that require advanced understanding of fitness in cancer.Longitudinal xenoengraftment interrogated via next generation single cell sequencing (SCS) hasenabled more accurate, quantitative measurements of tumours as they evolve. This process gen-erates timestamped samples comprising thousands of cells each measured at thousands of copynumber aberrations (CNA). Major analytical challenges introduced by this new datatype include(i) how to identify biologically meaningful groups of cells (i.e., clones) across multiple timepoints,and (ii) how to quantitatively reason about the underlying evolutionary forces acting on the clonesvia their observed dynamics.To address the first problem, we describe and provide supplementary tools for sitka, a scalableBayesian phylogenetic inference method in Chapter 2. It resolves the clonal structure of a heteroge-neous tumour cell population sampled over multiple timepoints by reconstructing the evolutionaryrelationship between single cells from their inferred CNA profiles. We then develop Lumberjack,a tree-cutting algorithm, and use it to assign cells to clones. We address the second problem inChapter 3 by developing fitClone, a Bayesian probabilistic framework that ascribes quantitativeselection coefficients to individual cancer clones and forecasts competitive clonal dynamics overtime. In Chapter 4, we exemplify the computational models introduced above on real-world datacollected from cancer cells over a multi-year period to verify two key hypotheses, that (i) clonaldynamics in a pre-treatment triple negative breast cancer (TNBC) tumour is quantifiably repro-ducible, and that (ii) the fitness landscape is reversed under early response to cisplatin treatment.Our results show that population genetic modelling of timeseries tumour measurements to predictiiiclonal evolution is tractable. Further study with timeseries modelling will provide insight into ther-apeutic strategies promoting early intervention, drug combinations and evolution-aware approachesto clinical management.ivLay SummaryCells accrue mutations that in time result in uncontrolled growth and in some cases the abilityto invade the surrounding tissue and seed metastasis. This thesis develops a statistical frameworkthat characterises the subpopulations that exit in an evolving tumour monitored over a period oftime at a single-cell resolution. This information is then used to generate a model that predictshow cancerous tumours evolve over time, with or without treatment. We expect that this researchprovides insight into therapeutic strategies promoting early intervention, drug combinations andevolution-aware approaches to clinical management of human cancers.vPrefaceA version of Chapter 2 is a manuscript under preparation and has a draft on the bioRxiv preprintrepository [1]. Professor Alexandre Bouchard-Coˆte´ developed and implemented the sitka model. Ideveloped the visualisation methods and performed the benchmarking experiments and generatedthe figures.Chapter 3 and Chapter 4 are based on a manuscript under review (bioRxiv preprint [2]). I havedeveloped and implemented the fitClone model and generated the figures in Chapters 2, 3, and 4.I have performed the phylogenetic and fitness analysis of all the datasets introduced in Chapters 2and 4. Farhia Kabeer has generated the xenograft biological systems. Dr. Marc Williams has donethe missegregation analysis. Professor Sohrab Shah and Samuel Aparicio have contributed to theoriginal writing of the manuscript on which these two chapters are based.viTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xivDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Tumour evolution and models of fitness . . . . . . . . . . . . . . . . . . . . . . . . . 21.1.1 Detection and measurement of evolutionary forces . . . . . . . . . . . . . . . 31.2 Measurement of clonal abundance over time . . . . . . . . . . . . . . . . . . . . . . 51.2.1 Leveraging single cell sequencing . . . . . . . . . . . . . . . . . . . . . . . . . 51.2.2 Clone identification via phylogenetic analysis . . . . . . . . . . . . . . . . . . 61.3 Research contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61.3.1 Chapter 2: sitka: benchmarking and visualisation of a scalable Bayesianphylogenetic inference method . . . . . . . . . . . . . . . . . . . . . . . . . . 71.3.2 Chapter 3: Modelling fitness in longitudinal data via an approximation tothe diffusion model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7vii1.3.3 Chapter 4: Application of the Wright-Fisher diffusion approximation to can-cer model systems interrogated at single cell level . . . . . . . . . . . . . . . 72 Benchmarking and Visualisation of a Scalable Bayesian Phylogenetic InferenceMethod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.2 The sitka model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.3 Workflow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.3.1 Filtering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.3.2 Binarising . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.3.3 Fixing jitter and selection of phylogenetic markers . . . . . . . . . . . . . . . 222.3.4 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222.3.5 Inference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 252.3.6 Summarising the posterior distribution . . . . . . . . . . . . . . . . . . . . . 272.3.7 Consensus tree and CNA heatmap visualisation . . . . . . . . . . . . . . . . 292.4 Overview of datasets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 332.5 Tree evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 352.5.1 Predictive test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 352.5.2 Comparison to baseline methods . . . . . . . . . . . . . . . . . . . . . . . . . 372.6 Identifying clones through phylogenetic analysis . . . . . . . . . . . . . . . . . . . . 412.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 423 Modelling Fitness in Longitudinal Data Via an Approximation to the Wright-Fisher Diffusion Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 453.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 453.1.1 fitClone: a Bayesian fitness model for timeseries data . . . . . . . . . . . . . 463.2 Wright-Fisher diffusions with selection . . . . . . . . . . . . . . . . . . . . . . . . . . 473.3 fitClone model specification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 483.4 Posterior inference under the fitClone model . . . . . . . . . . . . . . . . . . . . . . 503.4.1 Euler-Maruyama discretisation . . . . . . . . . . . . . . . . . . . . . . . . . . 503.4.2 Inference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51viii3.5 Conditional sampler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 523.6 Estimating the effective population size . . . . . . . . . . . . . . . . . . . . . . . . . 543.7 Summarising the posterior distribution . . . . . . . . . . . . . . . . . . . . . . . . . 553.8 Selecting the reference clone . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 563.9 The deterministic logistic growth model . . . . . . . . . . . . . . . . . . . . . . . . . 563.10 Simulation benchmarking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 573.11 Proof of principle application on TNBC-SA501 PDX . . . . . . . . . . . . . . . . . 593.11.1 Software and implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . 613.12 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 614 Application of the Wright-Fisher Diffusion Approximation to Cancer Model Sys-tems Interrogated at the Single Cell Level . . . . . . . . . . . . . . . . . . . . . . . 634.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 634.2 In vitro systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 654.2.1 Serial passaging of immortalized 184hTert diploid breast epithelial cell lines 654.2.2 Phylogenetic and fitness analysis . . . . . . . . . . . . . . . . . . . . . . . . . 664.2.3 Segmental aneuploidies correlate with positive selection in diploid p53 defi-cient cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 684.3 In vivo system in the untreated regime . . . . . . . . . . . . . . . . . . . . . . . . . 704.3.1 Serial passaging of TP53 mutant PDX tumours . . . . . . . . . . . . . . . . 704.3.2 Mixture experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 734.4 In vivo systems under cisplatin treatment . . . . . . . . . . . . . . . . . . . . . . . . 764.4.1 TNBC-SA609 PDX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 764.4.2 Fitness landscape reversal in early response to cisplatin . . . . . . . . . . . . 794.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 854.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 864.6.1 Limitations and future research . . . . . . . . . . . . . . . . . . . . . . . . . 865 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 885.1 Summary of contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 885.2 Future research directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90ix5.3 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93AppendixA.1 TNBC PDX experimental design diagrams. . . . . . . . . . . . . . . . . . . . . . . . 111A.2 Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113A.3 Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121A.4 Preprocessing of the CNA data for phylogenetic inference . . . . . . . . . . . . . . . 123A.5 Computing CNA rate over time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125xList of Tables2.1 The observation error matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 235.1 Repositories hosting the source codes of the developed softwares . . . . . . . . . . . 88A.2 Estimated CNA rate over time in timeseries datasets used in Chapter 2 . . . . . . . 121A.3 Parameters per dataset . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121A.4 Real-world data parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121A.5 Summary of real-world datasets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122xiList of Figures1.1 A diagram summarising the organisation of the thesis . . . . . . . . . . . . . . . . . 92.1 Description of the process involved in the construction of markers, the input to thesitka model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.2 Visualisation of a small tree and its corresponding marker matrix . . . . . . . . . . . 132.3 Perfect phylogeny and effects of overlapping CNA events. . . . . . . . . . . . . . . . 152.4 Workflow of sitka . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.5 Filtering the CNA data for tree inference. . . . . . . . . . . . . . . . . . . . . . . . . 182.6 Synchronous CNA changes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212.7 Pathological tree reconstruction under default observation prior . . . . . . . . . . . . 242.8 Marker-centric view for SA535 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 312.9 Phylogenetic tree and CNA profile heatmap for the SA535 dataset . . . . . . . . . . 342.10 Phylogenetic tree and CNA profile heatmap for the OV A dataset . . . . . . . . . . . 352.11 Phylogenetic tree and CNA profile heatmap for the SA501 dataset . . . . . . . . . . 352.12 Results over real-world datasets and benchmarking against baseline methods . . . . 403.1 A state space graphical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 493.2 Simulation studies for the fitClone model . . . . . . . . . . . . . . . . . . . . . . . . 583.3 The fitness analysis for clones derived from bulk WG . . . . . . . . . . . . . . . . . . 614.1 Schematic overview of experimental design for quantitatively modelling clone-specificfitness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 654.2 Impact of p53 mutation on modes of selection on 184hTert cells . . . . . . . . . . . . 694.3 Structural variant and mutation rates of 184hTert cells . . . . . . . . . . . . . . . . . 704.4 Comparison of fitness landscapes of breast cancer PDX models. . . . . . . . . . . . . 72xii4.5 Mixture experiments on TNBC-SA609 PDX . . . . . . . . . . . . . . . . . . . . . . . 754.6 Impact of pharmacologic perturbation with cisplatin on fitness landscapes in TNBC-SA609 PDX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 784.7 TNBC-SA535 PDX clonal dynamics with and without treatment . . . . . . . . . . . 814.8 TNBC-SA1035 PDX clonal dynamics with and without treatment . . . . . . . . . . 824.9 Fitness landscape reversal in early cisplatin treatment in TNBC PDX models . . . . 84A.1 The experimental design diagram for TNBC-SA609 PDX . . . . . . . . . . . . . . . 111A.2 The experimental design diagram for TNBC-SA535 PDX . . . . . . . . . . . . . . . 112A.3 The experimental design diagram for TNBC-SA1035 PDX . . . . . . . . . . . . . . . 113A.4 The distribution of mismatch rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113A.5 Credible interval for the selection coefficients of all clones across all datasets . . . . . 114A.6 Posterior order matrix for selective coefficients of all clones in all datasets . . . . . . 115A.7 Overview of experimental design and PDX growth curves . . . . . . . . . . . . . . . 116A.8 Single nucleotide variant (SNV) and copy number profiles of the in-vitro and in-vivosystems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117A.9 TNBC-SA609 PDX clonal dynamics with and without treatment . . . . . . . . . . . 118A.10 Structural variant and mutation rates of PDX lines . . . . . . . . . . . . . . . . . . . 119A.11 Fitness validation with tumour mixing and drug perturbation . . . . . . . . . . . . . 120A.12 An example of replicating cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124xiiiAcknowledgementsI would like to thank my supervisors, Dr. Sohrab Shah and Dr. Alexandre Bouchard-Coˆte´ for theirunyielding support both in my academic and personal life. Thank you to my supervisory committee,Dr. Samuel Aparicio and Dr. Sarah Otto. I learned a lot in terms of biological and evolutionaryfundamentals of the problems investigated in this work. I consider myself utterly fortunate to havehad a supervisory committee that were dedicated and engaged in my project.I would like to thank members of the Shah and Aparicio labs, specially Farhia Kabeer, KieranCampbell, Andrew McPherson, and Mirela Andronescu for their help and stimulating conversa-tions.xivDedicationTo my family, my gratitude for your unconditional support. To Monique, thank you for believingin me. Without you I could not have done this.xvChapter 1IntroductionCellular fitness underpins the tissue population dynamics of cancer progression and treatmentresponse. Yet, quantifying fitness in heterogeneous cell populations and identifying causal mecha-nisms shaping fitness landscapes remain open problems, impeding progress in developing effectiveand durable therapeutic strategies. In particular, quantitative fitness modelling of cancer cells hasnumerous and diverse implications: attributing clonal dynamics to drift or selection, identifyingthe determinants of clonal expansion, enabling causal inference, and forecasting growth trajecto-ries. Why and how drug resistance evolves is among the key unresolved areas of investigation thatrequire advanced understanding of fitness in cancer. For example, drug resistance mechanisms arecommonly attributed to phenotypic plasticity encoded via epigenetic changes [3, 4] or evolutionaryselection of pre-existing genomic clones [5]. However, the relative contribution of these processeswhen studied in tandem is poorly understood and requires integrated genome-transcriptome inves-tigation. Moreover, how changes in genomic architecture brought about by copy number alterations(CNA) drive tumour progression remains understudied [6]. Genetic editing of cancer related genesand treatment with pharmacological drugs are among possible ways to induce fitness changes. Theeffects of these perturbations can be monitored to profile tumour growth progression and othertherapeutic responses. We contend that quantitatively ascribing fitness values to clonal dynamicsover long-range timeseries in the context of such perturbations would provide higher order insightsin drug resistance than current models allow.Previous work has established models of fitness through interpreting allelic measurements ofsingle snapshots [7–11] from bulk sequencing over large patient cohorts [12], timeseries monitor-ing of cell free DNA [13], multiregion sequencing [8, 14–17] and estimating fitness landscapes ofclonal haematopoiesis [18]. However, the cancer field has generally lacked serial measurementsfrom patient derived tissues to directly observe cancer evolution over realistic timescales. This has1impeded a thorough understanding of factors driving selection, achieved in other biological systemsthrough studying granular timeseries with population genetic modelling [19]. The majority of workin cancer has focused on bulk tumour sequencing, where it is computationally difficult or sometimesimpossible to accurately resolve the subclonal composition of tumours. Single cell genome mea-surements to scalably define clonal populations in cancer over thousands of cells have only recentlyemerged [20, 21], enabling identification of rare populations, precise tracking of clones and robustclone-specific measurements suitable for population genetic modelling.In this dissertation we present a body of work that aims to (i) identify clones, (ii) measure theirabundances over time and (iii) infer evolutionary fitness parameters for populations measured at thesingle cell resolution to help investigate how human cancer cells evolve at the copy number level, andhow establishing baseline fitness measures helps to interpret selection under drug administration.1.1 Tumour evolution and models of fitnessCells accrue mutations that in time result in uncontrolled growth and in some cases the ability toinvade the surrounding tissue and seed metastasis. The stochastic nature of mutation accumulationand subsequent selection dictates that the tumour growth process is, to some degree, random. Thisin turn may result in treatment-resistant subpopulations and ultimately cause relapse. Advances innext generation sequencing (NGS) combined with statistical deconvolution methods have enabledinvestigation of tumour subpopulation structures at a high resolution [22, 23], where clusters ofmutations existing at the same abundance can be identified. Many resistance mechanisms todrugs that occur reproducibly in patients, and thus are predictors of response, are known [24, 25].For instance, the emergence of T790M mutations in EGFR in patients with non-small-cell lungcancers after treatment with tyrosine kinase inhibitors (TKI) [26], secondary mutation in KRASin colorectal cancer patients post mono-clonal antibody (mAb) anti-EGFR therapy [27, 28], ESR1mutations induced acquired-resistance in ER+ breast cancer patients undergoing endocrine therapy(e.g., tamoxifen) [29, 30] and BRCA1 revertant mutations in platinum and PARPi treated ovariancancers [31]. These illustrate that targeted selective pressures induce selection of specific alleles.Despite this, the literature is incomplete and we do not yet know what contributes to relapse inmost cases. Little progress has been made in the capacity to predict the genetic makeup of a2tumour sample over time and under drug administration.1.1.1 Detection and measurement of evolutionary forcesOne way to quantify the strength of selective pressure at the gene level is by measuring the rateof non-synonymous to synonymous mutations (dN/dS also known as Ka/Ks) [32]. Assuming thatthe infinite sites model holds, the probability of a mutation is uniform across the genome. Wewould expect that in the absence of any selection, dN/dS ≈ 1, since no functional change shouldbe advantageous or deleterious to the cell. dS (Ks) could be thought of as representing the rateof background mutation. Accurately estimating dN/dS is challenging. Multiple biases have to becorrected for, including the differential rate of mutation based on the upstream and downstreamnucleotides of a site (i.e., context), structural changes, and in the absence of normal matchingtissue, polymorphic sites. Additionally, very deep targeted sequencing would be required to ensurea reasonable number of mutations are recovered [33, 34]. A recent method demonstrating theseissues uses a tri-nucleotide context-aware substitution model with 192 parameters to overcome someof the mentioned biases, but it is only able to achieve statistical power at the level of cohorts ofpatients and not individuals, limiting its clinical feasibility [12].Another approach is to assume neutral evolution in an exponentially growing population whenthe inverse allele frequency, 1/f , has a linear relationship with the expected number of mutationsM [7, 35]. According to this framework, deviation from linearity between 1/f and M is evidenceof the presence of selection. An extension of this work allows for quantifying selection magnitudesfor up to two subpopulations [8]. Frequency of mutations are detected using bulk WGS of a singletimepoint, the state of the tumour at the time of a biopsy. The method in [7] is predicated upon theobservation that branching processes [36], a family of stochastic models intensely used in populationgenetics, in the presence of some modes of selection, fails to generate data that show a linear 1/fand M relation.A promising study concerns timeseries monitoring of cell free DNA in 45 patients with colorectalcancer and wildtype RAS, treated with Cetuximab [13]. Plasma samples were collected monthlyuntil disease progression state was determined. One application was to predict time to diseaseprogression. The tumour in each patient was partitioned into sensitive and resistant subpopulations.The tumour burden (the total number of cancerous cells) was modelled as the sum of the exponential3decay and the exponential growth of the sensitive and resistant subpopulations respectively, eachparameterized by the initial population size and a decay (growth) rate. A key promise of thisstudy is the relative ease of obtaining plasma samples from patients as compared to solid biopsies.This allows more frequent sampling and forming of longitudinal timeseries. A drawback is thatthe sequencing of the plasma samples was limited (a panel of predefined 77 oncogenes and genomicregions), which in turn hindered the ability to detect all existing subpopulations. The results werealso focused on one type of cancer and the generalisability to other cancer types is unknown.An experimental system, the CRISPR-knockout screens [37, 38], allows the abundance of sub-populations (genotypes) of isogenic backgrounds to be observed, where each harbours specific mu-tations induced by CRISPR-Cas9. Including a set of non-targeting sgRNAs establishes a baselinefor neutral growth, i.e., the expected change in the abundance of putatively neutral genotypes. Insuch assays, depletion or expansion of a genotype hints at negative or positive selection respectively.An improved version of these screens incorporates unique molecular identifiers (UMIs) [39] into thesgRNAs allowing for quantitative interpretation of the effective population size of each expandedgenotype, unaffected by PCR duplicates [40].Inferring population genetics parameters with high confidence based on only one timepoint isdifficult and the associated uncertainty may render making meaningful conclusions implausible [41].For one, as shown in [41], models of evolution with selection can also produce a linear 1/f and Mrelationship. Birth, death and mutation rates, population size and number of generations are allunknown and comprise too many degrees of freedom to fit to one noisy observation of the endpointof a stochastic process. The signal for clones existing at similar prevalences is obscured in bulkWGS data and the frequency of alleles is further skewed by CNAs, which are difficult to resolveusing bulk WGS data [23, 42]. Analysis based on a single sample is particularly prone to samplingbias due to the spatial heterogeneity of the tumour and may result in loss of signal. Low sequencingdepth, characteristic of public libraries denies access to more recent events which further dilutes thesignal. These limitations combined suggest that bulk WGS at a single-timepoint may be insufficientto quantify fitness values with high confidence.41.2 Measurement of clonal abundance over timeLongitudinal studies, including immortalised cell lines and xenoengraftment experiments, allowtracking changes in the relative abundances of clusters of mutations [43], as well as Rx dosage ex-periments to induce resistance with subclinical dosing. These dynamics appear to be reproducible,raising hope that they may be predictable. Patient derived xenograft (PDX) systems are an effec-tive model to study timeseries of a human tumour cell population [44, 45]. By serially transplantinga patient derived tumour sample into highly immunodeficient mice, it is possible to continuouslymonitor how tumour composition may evolve in a patient over time. CRISPR-Cas9 has made di-rected evolution experiments possible in both cell lines and PDX systems [40, 46]. Specific genomiccontexts could be set up to study gene interactions in the presence of genetic perturbation inducedby CRISPR-Cas9.1.2.1 Leveraging single cell sequencingAdvances in SCS permits the investigation of tumour genetic composition with unprecedentedaccuracy and resolution. SCS obviates the need for complex computational deconvolution methodsand allows the direct observation of the genotypes of the extant subclones. This in turn allows theuse of models from the well established discipline of population genetics to be adapted to powerquantitative reasoning about the future trajectory, dominance, and depletion of subpopulations.Reproducible clonal dynamics from serially propagated PDX models with validated clonal genomesat single cell resolution has been established [43]. Motivated by this observation, we propose toinvestigate the hypothesis that clones that grow reproducibly have higher fitness.Real-time monitoring of evolution of cancerous cells at the single cells CNA resolution is nowfeasible [20, 21]. These advances in single cell whole genome sequencing enables population geneticsanalysis of cancer progression. It has been shown that single cell sequencing is essential for resolvinggenomically defined clones in high-grade serous ovarian cancer (HGSC) [17] and in patient derivedbreast cancer xenografts [43]. These studies employed multiplexed PCR approaches on amplifiedmaterial targeting small sets of mutations assayed with Fluidigm Access arrays. To enable wholegenome single cell sequencing, bias limitations in genomic depth [47] and breadth [48] have beenmitigated. A direct tagmentation approach [21] (direct library prep (DLP)) implemented in a micro-5lithography microfluidic device obviates pre-amplification and results in a superior representationof the copy number architecture of single cancer cells to below 1Mb resolution [20].Quantitative attributes of fitness are central to developing predictive models of cancer evolution.Such models currently elude the field due to a lack of appropriate timeseries measurements andgenerative analytic methods. The goal is to predict unobserved clonal trajectories over time byinferring the evolutionary dynamics underlying disease progression.1.2.2 Clone identification via phylogenetic analysisA key step in determining the quantitative fitness attributes of cancer cells is to resolve the sub-population structure. Cancer cells are evolutionarily related and phylogenetic reconstruction aimsat recovering this relationship. A fundamental assumption is that cells that descend from the sameancestor harbour similar mutations. While single cell whole genome sequencing (scWGS) has anadvantage over bulk WGS in that it obviates the need to computationally determine which mu-tations co-occur, the evolutionary history of the cells is not directly observed. Sorting cells intogroups enables us to track the change in their abundance over time. A principled way to groupcells is according to their phylogenetic relationships. Many phylogenetic inference methods aretailored for using point mutations as input and assume a small number of leaf nodes [49]. However,the emerging single cell data produce up to thousands of single cell genomes and are suitable fordetermining CNAs [20, 21].1.3 Research contributionsThis dissertation is organised as follows. In Chapter 2, we review and benchmark sitka, a scalableBayesian phylogenetic tree inference method and introduce Lumberjack, a tree-cutting algorithmto identify clones from sitka inferred trees. In Chapter 3, we introduce and develop the fitClonemodel, a Bayesian fitness inference framework for timeseries data and establish its performance viasimulation studies. Finally, in Chapter 4, we apply methods developed in the previous chapters ontimeseries data from in vitro and in vivo model systems with and without treatment. Figure 1.1shows a diagrammatic representation of this dissertation. Below we give a summary of the chaptersto come.61.3.1 Chapter 2: sitka: benchmarking and visualisation of a scalable Bayesianphylogenetic inference methodA resolved clonal structure of a heterogeneous tumour cell population sampled over multiple time-points is a key input to the fitClone model. In many modern single cell DNA sequencing scDNAseqplatforms including DLP+ used in this work, the evolutionary relationship between single cells isnot directly observed. sitka is a scaleable method that approximates the phylogeny in a Bayesianframework. In this chapter we first review the sitka probabilistic model and its inference pro-cedure. We then describe its workflow, starting from preprocessing its inputs to visualising itsoutputs. Next we exemplify sitka on three real-world cancer datasets, namely two from TNBC-PDX timeseries and one multi-region pre and post treatment HGSC sample. We review someexisting baseline methods and compare them to sitka using an accuracy test that we derive. Fi-nally we introduce Lumberjack, a method for cutting the phylogenetic tree to assign cells to clonesdetermined by their CNA profiles.1.3.2 Chapter 3: Modelling fitness in longitudinal data via an approximationto the diffusion modelIn this chapter we develop a probabilistic approach that is a mechanistic model of tumour growthbased on the Wright-Fisher diffusion approximation. We introduce the fitClone probabilisticframework, its observation model and inference procedure. We then describe an extension, theconditional sampler, that would allow inference for a larger number of clones. We review parameterestimation and key posterior summarisation procedures, specifically the posterior ordering matrixthat can be used for pairwise comparison of clones and in a principled way establish whether aclone has higher fitness than another. We conclude by exemplifying the method over simulationstudies and a real-world dataset.1.3.3 Chapter 4: Application of the Wright-Fisher diffusion approximation tocancer model systems interrogated at single cell levelIn this chapter we exemplify the computational framework that we have developed in previouschapters on real-world dataset to make inferences about how human cancer cells evolve at the7copy number level. We address two key hypotheses, that (i) clonal dynamics in a pre-treatmenttriple negative breast cancer (TNBC) tumour is quantifiably reproducible, and that (ii) the fitnesslandscape is reversed under early response to cisplatin treatment. For all datasets, we use sitkaand Lumberjack to infer phylogenetic trees and assign cells to clones, and then use fitClone toascribe quantitative selection coefficients to individual cancer clones.As a proof of principle, we first look at the dynamics of three related cell lines. We thenexamine two independent TNBC PDX samples and in one, TNBC-SA609, we validate our selectioncoefficient predictions in-vivo via a competition experimental design that addresses hypothesis (i).To tackle hypothesis (ii) we look at three independent TNBC PDX systems, with and withouttreatment. We then explore in more detail the dynamics of the TNBC-SA609 system. We alsolook at a drug holiday experimental design, in which treatment is stopped after a few timepoints,that suggests that withholding treatment may reset the clonal composition of the tumour to apre-treatment state.8Describing sitka single-cell Bayesian phylogenetics Preprocessing Visualisation BenchmarkingChapter 2Clone identification viaPhylogeneticsfitClone: Modelling fitness in timeseriesProbabilistic modelInference engine and extensionSimulation & real-world data benchmarkingChapter 3Ascribing quantitative selection coefficients to individual clonesLumberjack: Cutting phylogenetic treesThe fitness landscape is reversed under early response to cisplatin treatmentInterpret selection under drug administrationFitness landscape reversal TreatedSingle cell fitness landscapesIn-vitroExemplify the analytical frameworkIn-vivoChapter 4Application tomodel systemsCompetition mixture experimentsUntreatedClonal dynamics in a pre-treatment TNBC tumour is quantifiably reproducibleVerification of selection coefficientsFigure 1.1: A diagram summarising the organisation of the thesis.9Chapter 2Benchmarking and Visualisation of aScalable Bayesian PhylogeneticInference Method2.1 IntroductionA main challenge in investigating cancer evolution is the need to resolve the subpopulation struc-ture of a heterogeneous tumour sample. Advances in next generation SCS as well as longitudi-nal xenoengraftment has enabled more accurate, quantitative measurements of tumours as theyevolve. Phylogenetic reconstruction is central to identifying clones over multiple sampled time-points. Moreover, the ability to sequence tens of thousands of single genomes at high resolutionper experiment [20] calls for tailored phylogenetic models with scalable inference algorithms.Single cell cancer phylogenetics is an evolving field. Multiple approaches, spanning differentstudy designs and data sources are reviewed in [50]. Many phylogenetic inference methods aretailored for using point mutations as input and assume a small number (on the order of 10s) ofleaf nodes [49, 51–53]. However, the emerging single cell data produce up to thousands of singlecell genomes and are suitable for determining CNAs [20, 21]. Distance based and agglomerativeclustering methods such as neighbour joining are scalable and are used to elucidate hierarchicalstructures over cells [54, 55]. While useful heuristics, these methods are statistically sub-optimalrelative to likelihood based methods [56].We describe sitka, a phylogenetic model and associated Bayesian inference procedure whichexploits the specifics of scWGS data. sitka enjoys a number of attractive characteristics, including(i) a novel phylogenetic encoding of CN data providing a statistical-computational trade-off by10simplifying the site dependencies induced by rearrangements while still forming a sound foundationto phylogenetic inference, and (ii) an innovative phylogenetic tree exploration move which makesthe cost of MCMC iterations bounded by O(|C|+|L|), where |C| is the number of cells and |L| is thenumber of loci. In contrast, existing off-the-shelf likelihood-based methods incur an iteration cost ofO(|C| |L|), and (iii) the novel move considers an exponential number of neighbouring trees whereasoff-the-shelf moves consider a polynomial size set of neighbours. A key piece is to add visualisationtools to the sitka toolbox and this is a main contribution of this chapter. Visualisation techniqueshelp explore, summarise and communicate the tree inference results.In this chapter, we will first give a high level description the sitka model (Section 2.2) andreview its analysis workflow (Section 2.3). We will present the quality control (QC) of singlecells, different encodings of the input space (preprocessing), and then formally describe the model(Section 2.3.4) and an MCMC based inference procedure (Section 2.3.5). We will introduce twonovel visualisation views to summarise the tree inference results (Section 2.3.7). We will comparesitka with other tree-inference methods that scale to scWGS datasets over three real-world datasets(Sections 2.4 and 2.5). Finally we introduce Lumberjack, an algorithm for cutting a phylogenetictree to determine clonal subpopulations (Section 2.6).2.2 The sitka modelsitka is based on lossy transformation of single cell copy number matrices retaining only presenceor absence of changes in copy number profiles. This transformation turns a complex evolutionaryprocess (integer-valued copy numbers, prone to a high degree of homoplasy and dense dependencestructure across sites) into a simpler one which can be approximated by a probabilistic version of aperfect phylogeny (see Figure 2.1). We leverage the special structure created by the change pointtransformation to build a special purpose MCMC kernel which has better computational scalabilityper move compared to classical phylogenetic kernels.The input data for sitka can be visualised in a colour-coded matrix exemplified in Figure 2.1-a.Each row in the matrix corresponds to an individual cell that has been sequenced in a single-cellplatform. Each column in the matrix is a locus that is represented by a bin, a contiguous setof genomic positions. We assume that the integer copy number of each bin has been estimated11CNA matrixcell 1cell 2cell 3cell 4cell 5bin 1bin 2bin 3bin 4a b{A chromosome062817 93 410511+Copy number state062817 93 410511+Copy number statecell 1cell 2cell 3cell 4cell 5marker 2marker 3marker 1marker matrixcell 1cell 2cell 3cell 4cell 5bin 1bin 2bin 3bin 4marker 2marker 3marker 1CNA interlaced with markerscMarkerPresentAbsentMarkerPresentAbsentFigure 2.1: Description of the process involved in construction of markers, the input to the sitkamodel. A bin is a contiguous set of genomic positions. Each pair of consecutive bins (e.g., bins 1and 2 in (a) is associated with a marker (e.g., marker 1) that measures for each individual cell,whether there is a difference between the CNA states of the two bins. (a) The observed CNA matrixfor a subset of bins on a chromosome. The rows are sequenced single cells, and the columns arebins. The CN states are colour-coded. (b) The three markers shown are associated with the fourbins. Each marker records the presence (black) or absence (white) of a CN state change betweena pair of consecutive bins. Note that in the CNA matrix, there is a CN change at row 3 from bin1 to bin 2 (CN state 3 to 6). This is reflected in the marker matrix, at row 3 of marker 1 witha black square. There are no changes between bins 2 and 3 across any rows in the CNA matrix.This is reflected in marker 2 comprising all white squares. (c) For visualisation purposes, the CNAmatrix can be interlaced with the marker matrix to more clearly show where the CNA changesoccur. Each column of the marker matrix is inserted between the associated pair of columns in theCNA matrix. The resulting matrix is an example of an augmented view that combines data fromtwo or more sources (here the CNA matrix and the marker matrix). In an augmented view, we callcolumns from each source a channel.(i.e., called) as a preprocessing step, e.g., using a hidden Markov model tool such as in [21]. InFigure 2.1-a the copy number state is encoded by the colour of each entry in the matrix.The output of sitka includes two types of trees. The first is the tree used for MCMC samplingin the inference procedure (type I), and the second is the tree used in visualisation (type II) andcan be derived from the first type. Briefly, the first type is a directed tree that spans all cells andCN change points (markers) under study in addition to a virtual root node. If a cell vertex isattached to a marker vertex, we hypothesize that the cell harbours only markers in the shortestpath between its parent and the root node in the tree. A cell vertex, attached to the root node,is hypothesized to harbour none of the markers. If a marker vertex is a descendent of another12marker vertex, we hypothesize that the trait that the former indexes, emerged in a cell that was adescendant of the cell that the latter indexes. Note that in this type of tree, some marker verticescould be leaf nodes. Such markers are hypothesised to violate the assumptions of the model andwere not consistent with any cells. All cell vertices will be leaf nodes. We describe the type I treemore formally in Section 2.3.4 when we introduce priors on trees.Type II tree is a transformation of type I tree. We remove from the type I tree all marker nodesthat are leaf nodes, i.e., markers that are not present in any cells. We also collapse into a singlenode, the list of connected marker nodes that have exactly one descendent (i.e., chains). Figure 2.2shows a small type I tree, its transformation to a type II tree and the respective marker matrix.We visualise the input matrix and the estimated tree simultaneously by sorting the individual cells(rows of the matrix) in such a way that they line up with the position of the corresponding leavesof the tree.marker 2marker 3marker 1Latent marker matrix Figure 4: Visualization of a small tree t 2 T . The right side shows the matrix of binary latentvalues xc,l where rows are cells and columns are loci. The left side shows a summary of the perfectphylogeny. To remove clutter, we remove from the tree (but not the matrix) all loci where xc,l = 0for all cells. We also collapse chains (list of edges with exactly one descendant) into a single edge.Note that all clade containing more than one descendant cells is necessarily supported by one ormore traits l having value xc,l = 1 if and only if c is a descendant of the clade. Clades at the leaves,which have size one may or may not be supported.Sampling a mutation tree: let Vm = L [ {v⇤} denote a vertex set composed of one vertex foreach of the |L| loci plus one artificial root node v⇤. The artificial root node induces an implicitnotion of direction on the edges, viewing them as pointing away from v⇤. Let T m denote theset of trees tm spanning Vm. The interpretation of tm is as follows: there is a directed pathfrom vertex/locus l to l0 in tm if and only if the trait indexed by l is hypothesized to haveemerged in a cell which is ancestral to the cell in which l0 emerged.Sampling cell assignments: assign each cell to a vertex in tm. The interpretation of assigningcell c to locus l is that among the traits under study, c is hypothesized to possess only the traitsvisited by the shorted path from v⇤ to l in tm. If a cell c is assigned to v⇤, the interpretationis that c is hypothesized to possess none of the traits under study.Both steps can be viewed as graphs (for the second step, it is a bipartite graph with onecomponent being the set of loci, and the other, the set of cells). It is convenient to summarize bothprocesses at the union of the two graphs, which is also a tree, this time on |L| + |C| + 1 vertices,V = L [ C [ {v⇤}. Let us denote the set of trees obtained by this two step process by T . Again,the vertex v⇤ induces an implicit direction to edges in t 2 T .Given t 2 T , the matrix x is a deterministic function obtained by setting xc,l = 1 if vertex c isa descendant of vertex l in t, and zero otherwise. We denote this deterministic function by x(t).We show a small example in Figure 4.The tree structure encoded by T is unidentifiable: for example, if a tree contains a chain of traits,then permutations of the traits yield the same matrix x(t). Since we are taking a Bayesian approachthis is not problematic. As we describe later, tree summaries we build from the posterior distributioncollapse the non-identifiable parts. Moreover, the specific choice we make for constructing the setT help us obtain simple and ecient sampling algorithms.To complete the construction of the prior on trees, we need to assign probabilities to elements9MarkerPresentAbsentcell 1cell 2cell 3cell 4cell 5Typ I t eerootcell 1cell 2cell 3cell 4cell 5Type II tr ea b cCell Marker RootVertex typeFigure 2.2: Visualisation of a small type I tree t (a), its transformation into a type II tree (b),and the corresponding marker matrix x = (xc,l) (c). Given a tree t, the latent marker matrix xis a deterministic function x = x(t). We compute x : t → {0, 1}C×L by setting xc,l = 1 if thesingle-cell c is a descendent of the marker node l in tree t, and otherwise xc,l = 0. Note thatthe clade comprising single-cells 3 and 4 has support in both markers 1 and 3. For clarity, we donot visualise type I trees, but plot their transformation, i.e., type II trees as follows. We removefrom the type I tree ll marker nodes that have xc,l = 0 for all single-cells c. List of connectededges that have exactly one descendent (i.e., chains) are also collapsed into a single edge, e.g., theedge corresponding to markers 2 and 3 are collapsed into one edge (since marker 2 has only onedescendent, namely single-cell 2).sitka is predicated on the perfect phylogeny assumption where each phylogenetic trait arisesexactly once on a rooted tree topology and that all individual cells descending from that positionon the tree will inherit that trait. The perfect phylogeny model allows us to use CNA change13points, i.e., markers, instead of integer CNA states as phylogenetic traits. When using changepoints as phylogenetic traits, non-overlapping CNA events do not break the perfect phylogenyassumption. Figure 2.3 shows examples of overlapping CNA events and their effect on markers.The two scenarios that can lead to the violation of the perfect phylogeny assumption are (i) whena CNA gain event is followed by an overlapping loss event, or (ii) when a loss event is followed byan overlapping loss event, and the second event removes either end-point of the first event. Notehow for a violation to occur, the second overlapping event has to happen on the same copy as thefirst event.Imposing a perfect phylogeny on the observed change points is restrictive, as violations of theassumptions (e.g., due to homoplasy) or an error in observing a single trait can lead to considerablechange to the inferred tree. To improve the robustness of the perfect phylogeny construction, we usean observation model described in Section 2.3.4. A central modelling assumption is that the perfectphylogeny holds for the latent states (denoted by xc,l = xc,l(t) in Section 2.3.4), but not for theobserved markers (yc,l in Section 2.3.4). Using this formulation, incorrect change point estimates (asa result of noisy CNAs calls) as well as markers inconsistent with the perfect phylogeny assumptioncan be accounted for as noise in the observations.2.3 WorkflowWe begin by describing the preprocessing and lossy transformation of the input CNA states andthen briefly describe the inference procedure. Figure 2.4 shows the workflow of the phylogeneticinference in sitka.14MarkerPresentAbsentbin 1bin 2bin 3bin 4marker 2marker 3marker 1Cellancestral state 2 2 2 2bin 5marker 4CellCNA gain event 2 3 3 3CellCNA gain event 2 32224 3aPerfect phylogeny maintained by two overlapping CN gains.c MarkerPresentAbsentbin 1bin 2bin 3bin 4marker 2marker 3marker 1Cellancestral state 2 2 22 2bin 5marker 4CellCNA gain event 2 2 23 2CellCNA loss event 2 1 21 1Perfect phylogeny violated by a CNA gain followed by a CNA lossbMarkerPresentAbsentbin 1bin 2bin 3bin 4marker 2marker 3marker 1Cellancestral state 2 2 22 2bin 5marker 4CellCNA loss event 2 2 21 2CellCNA loss event 2 2 21 1Perfect phylogeny violated by a CNA loss followed by a CNA lossViolation Violation062817 93 410511+Copy number state062817 93 410511+Copy number state062817 93 410511+Copy number stateFigure 2.3: The effects of overlapping CNA events on the perfect phylogeny assumption. A segmentof a chromosome with five consecutive bins and their four corresponding markers are shown. Eachpanel follows the CN states interlaced with markers for a cell at the ancestral state (top), after aCNA event (middle), and after a second overlapping CNA event (bottom). The numbers in theCNA squares show the integer CN state (e.g., the ancestral state has two copies of the 5-bins longsegment). (a) Two overlapping CNA gains maintain the perfect phylogeny assumption. By theinfinite site argument, it is unlikely for the end-points of the two gain events to exactly match.The same argument holds for a CNA loss followed by a CNA gain event. Note that in these cases,once a change point is acquired, it is not lost. (b) If a loss event is followed by another loss eventin which either end-points of the first event is removed, the perfect phylogeny assumption will beviolated (e.g., marker 3 is lost after the second loss event). Note that a violation does not occur ifthe loss events hit different copies of a segment. (c) Similarly, if a gain event is followed by a lossevent, only if the latter erases the end-points of the former is the perfect phylogeny violated. Notehow marker 2 and marker 3 are lost after the second CNA event.15| | | | | | | | | | || | | || | | | | | | || | | | | | | | | | || | | || | | | | | | |Copy numbers Observed markers Observed markersJitter fixed| | | | | | | | || | || | | | | | | || | | | | | | | | | || | | || | | | | | | || | | | | | | | | | || | | || | | | | | | || | | | | | | | | | || | | || | | | | | | |Selection of phylogenetic markers| | | | | | | | | | || | | || | | | | | | || | | | | | | | | | | | | || | | | | | |MCMCchr15_2950Iter = 100chr15_2950Iter = 101●●●R54C68R05C49R07C09chr17_2100chr15_2950chr1_4900iter = 101●●●R05C49R07C09R54C68chr12_1600chr15_2950iter = 100chr17_2100chr15_2950iter = 10Iter = 100 Iter = 101Iter = 1000…| | | | | | | | | | || | | || | | | | | | || | | | | | | | | | || | | || | | | | | | |AveragingEdge insertion| ||| |||| | || || | || | | || | | || | || | || || | | || | | || | ||| || || | || ||| | || | || | || || | ||| | | || || || ||| | ||| || | ||| | | |||| | || ||| || ||| || || ||| | || ||| || | | || || ||| | || |||| | | || | ||||| ||| | || | |||| | || | | || || || | | || | || | | || ||| ||| || ||| | | || | | || | ||| || | | || || | || | ||| | ||| || |||| || | ||| || ||| | | || ||| ||| | | |||| | || || | ||| | ||| | ||| || || |||| || | | ||| ||| | | | || || ||Consensus tree & sorted CNA heatmapa bd e fg h iCellsChromosomesCellsMarkers| | | | | | | | | | || | | || | | | | | | || | | | | | | | | | || | | || | | | | | | |cCellsMarkersCellsMarkersCellsMarkersCellsChromosomesCellsMarkerPresentAbsentMarkerPresentAbsentMarkerPresentAbsentSitka posteriorPresent0.0 1.0Inferred markerPresentAbsentSitka posteriorPresent0.0 1.0Observed markerPresentAbsent| || |22 XCN bins{ {Marker tuples (inf., post., obs.)062817 93 410511+Copy number state062817 93 410511+Copy number state062817 93 410511+Copy number state16Figure 2.4 (previous page): Workflow of sitka. (a) The input is the CNA data from a heterogeneoussingle-cell population. The rows of the CNA matrix are permuted at random. (b) CNA changepoints, i.e., markers are obtained by applying a lossy binary transformation to the CNA matrix(Section 2.3.2 and Figure 2.1). Note that each single-cell is now represented by the presence orabsence of CN changes between consecutive bins. (c) A shared change point is sometimes shiftedby a few bins across different cells. We use a heuristic (Section 2.3.3) to correct for this markermisalignment. Note how the columns in the inset in panel-c are less noisy than their counterpartin panel-b. (d) As a computational trade-off, only a subset of markers present in at least 5% ofcells are chosen for phylogenetic inference. (e) An MCMC algorithm explores the tree space. (f)An example of an edge insertion MCMC move. (g) The posterior probability of latent state xc,lis computed by averaging over all MCMC trees, summing over the number of samples in which asingle-cell c was a descendant of a marker node l (Section 2.3.6). (h) The rows in the CNA matrixin panel-a are sorted according to the inferred consensus (type II) tree, shown on the left of thematrix. Note that the block structure of the CNs are clearly visible as single-cells with similar CNprofiles are grouped together. (i) The inset shows tuples of marker columns, namely inf. (inferredmarkers, i.e., latent state xc,l), post. (posterior probability of the latent state xc,l), and obs.(observed markers), interlaced with the CN columns (similar to Figure 2.1). The results are fromthe SA535 dataset, a triple negative breast cancer patient derived xenograft sample (Section 2.4).172.3.1 FilteringThe raw data contains cells that are either contaminated (e.g., contains biological material frommice) or have undesired sequencing artefacts. These include cells that were captured for DNAsequencing when undergoing mitosis. Since the stka model does not account for such phenomena,the filtering is an important step. Figure 2.5 shows the steps taken from pulling the raw data to theCNA integer matrix ready for sitka transformation. See Section A.4 for detailed descriptions ofeach step in the filtering process. Briefly, we remove control cells, cells with highly-noisy CN calls,and cells that have very few mapped reads. We also remove copy number bins that lie in difficultto sequence regions of the genome (bins with low-mapability). Finally, we drop cells that, basedon their CNA profile, are suspected to be cycling cells that were not detected in previous steps.cn.csvCopy number and cell meta-datacn_bin_filtered.csvdrop low-mapability bins cn_bin_cell_filtered.csvdrop low-quality, contaminated, and cycling cellscn_bin_cell_filtered_no_jump.csvdrop suspect cycling cellsFigure 2.5: Filtering the CNA data for tree inference.2.3.2 BinarisingTo obtain the C × LMarkers phylogenetic markers matrix y that comprises the input to the sitkamodel, we apply a lossy transformation to the C × LBins CNA matrix a that involves computing18the change in copy number state between two consecutive bins. Below we introduce two differentways to compute the binarisation which differ on whether any artificial bins are included or not.Figure 2.1 shows a small CNA matrix and its corresponding transformation into the marker matrix.For brevity, in what follows we assume that only one chromosome is used, so that LBins = L andLMarkers = LBins − 1. In practice, we use all available chromosomes, and LMarkers = LBins − NChrwhere NChr denotes the total number of chromosomes used.BasicIn the basic binarisation process, the markers are computed as follows. Given a filtered cell-by-locusmatrix a (cn bin cell filtered no jump.csv), we sort bins by (chr-name, start-position,end-position). Then in each chromosome, we compute markers as the binarised difference betweenconsecutive bins, and drop duplicated markers, i.e., keep one marker in a set of markers withidentical binary patterns across cells. In other words, y = (yc,l′) and l′ ∈ {1, . . . , L− 1}, andyc,l′ := 1(∣∣ac,l′ − ac,l′+1∣∣ > 0) (2.1)where 1(x) is the indicator function.SyncSome datasets exhibit synchronous copy number changes in which at least two distinct cell subpop-ulations have a copy number change at a bin across a subset of cells, where either the direction orthe magnitude of the change is different in at least one cell. That is, there exists i, c, c′ such that:(yc,i = yc′,i = 1) ∧ (ac,i − ac,i+1 6= ac′,i − ac′,i+1) (2.2)for i ∈ {1, . . . , L− 1} and c, c′ ∈ {1, . . . , C} and a and y are the CNA and marker matrices respec-tively. Figure 2.6 shows examples of synchronous copy number changes. The basic binarisationmethod will be unable to distinguish such copy number induced subpopulations. The Sync proce-dure aims to alleviate this by inserting the following additional marker columns in matrix y (theresult of the basic encoding above):191. Find columns b := a1:C,i that satisfy Equation (2.2).2. For each b found above, add a virtual column b′ with its CN state set to the most frequentCN state in b.3. Add to the matrix y a column e whose elements are computed as the binarised differencebetween b and b′, that is, for c ∈ {1, . . . , C} set ec,1 := 1(∣∣∣bc,1 − b′c,1∣∣∣ > 0).20cell 2 has an identical set of change points to cell 3 in this chromosome, while CN states of the two cells are different over all bins. This can happen as a result of a whole chromosome duplication event.Pathological case 1. missing sync. CNAabin 1bin 2bin 3bin 4marker 2marker 3marker 1cell 1cell 2cell 3cell 4cell 5b{One chromosomePathological case 2: missing whole chr. eventMarkerPresentAbsent062817 93 410511+Copy number state062817 93 410511+Copy number stateCN states interlaced with markerscell 1cell 2cell 3cell 4cell 5bin 1bin 2bin 3bin 4marker 2marker 3marker 1MarkerPresentAbsent{One chromosomecell 3 and cell 4 have a copy number change frombin 1 to bin 2 (shown in black at marker 1), but the magnitude of the changesare different. CN states interlaced with markersFigure 2.6: Synchronous CNA changes. If two distinct cell subpopulations have a copy numberchange at a bin across a subset of cells, where either the direction or the magnitude of the changeis different in at least one cell (a synchronous CNA change), the basic binarisation method willbe unable to distinguish the induced subpopulations. Black and white represent the presence orabsence of a CN change (marker) respectfully. Panels (a) and (b) show cases where synchronousor whole chromosome events are missed.212.3.3 Fixing jitter and selection of phylogenetic markersThe copy numbers available to us in this work are estimated independently for each cell. Thisis one reason why the start position (bin) of the same CN change event may be slightly differentacross cells, generating some jitter. We address this by enumerating each change point column andmoving the jitter back to the most frequent bin. Concretely,Algorithm 1 JitterFix1: procedure jitter-fix(y, k)2: column-queue ← OrderByDensityDecreasing(y)3: columns-visited ← {}4: for column-index c in column-queue do5: neighbours ← neighbours (c, y, k)6: for column-index n in neighbours do7: if n /∈ columns-visited then8: y1:C,c ← y1:C,c ∨ y1:C,n9: y1:C,n ← 010: columns-visited ← columns-visited ∪ n11: return yThe function neighbours (i, y, k) returns 2k columns (if they exist), half, immediately to theleft, and half, immediately to the right of column i in matrix y. We use k = 2 based on visualinspection. An example of the result of the jitter correction heuristic is shown in Figure 2.4 panel c.As a computational trade-off, only a subset of markers present in at least a minimum number ofcells are chosen for phylogenetic inference. That is, we removed columns l in y with relative density∑c∈C yc,l/|C| less than a threshold, set to 5%. Larger values of this threshold may lead to lessresolved clades in the inferred tree.2.3.4 ModelThe sitka model starts with the perfect phylogeny assumption for the latent variables xc,l butallows deviation from it by allowing noisy observations yc,l. In a perfect phylogeny model, eachphylogenetic trait arises only once on the rooted tree topology, and all cells descending from thatposition will inherit that trait. Recall that we assume the latent states matrix is a deterministicfunction of the tree t, i.e., x = x(t) (see Figure 2.2).Let C and L denote the disjoint sets of cells and loci respectively. Assuming that a∗ is thetrue CNA state matrix, then informally, we can think of the trait xc,l for cell c ∈ C at markerl ∈ L as in Equation (2.1), with a replaced with a∗. But a∗ is not directly observed and therefore22neither is xc,l. We assume noisy estimates of CNA calls, ac,l are available. Then we define yc,l asin Equation (2.1). Informally, yc,l is a noisy version of xc,l. Therefore, we posit an observationprobability model p(y |x, θ), where x is the matrix x = (xc,l), y is the data or a data summary, andθ are model parameters. To model errors in copy number states, we introduce false positive andnegative rate parameters rFP ∈ (0, 1) and rFN ∈ (0, 1) respectively, and an error matrixerFP,rFN =1− rFP rFPrFN 1− rFN , (2.3)p(yc,l |xc,l, rFP, rFN)= erFP,rFNxc,l,yc,l, (2.4)from which we set:p(y |x, θ) =∏l∈L∏c∈Cp(yc,l |xc,l, rFPc,l (θ), rFNc,l (θ)),where the scalar erFP,rFNxc,l,yc,l takes values as in Table 2.1. Here we use a global parameterisation wherethe false positive and false negative functions rFPc,l (θ) and rFNc,l (θ) are shared across markers andcells. We considered a marker-specific parameterisation, but it incurs a higher computational costand results in a similar performance and that is why we focus here on a global parameterisation.xc,l yc,l pc,l0 0 1− rFP0 1 rFP1 0 rFN1 1 1 - rFNTable 2.1: The false positive and false negative errors in sitka’s observation model.In the global parameterisation, we have θ = (rFPglobal, rFNglobal). Using a uniform prior distributionfor both error rates can lead to pathological cases as shown in Figure 2.7. To avoid that, we usethe following:rFPglobal ∼ Uniform(0, rFP),rFNglobal ∼ Uniform(0, rFN),23We set rFP = 1/10 and rFN = 1/2 as defaults in our experiments. We use a larger rFN;although errors in approximating copy number states can lead to both types of errors, there aremore mechanisms that can cause a false negative, including a copy number gain, followed by anoverlapping copy number loss.✽✽ ✽✽✽✽ ✽True treelatentmarker 1latent marker 2cell 1 cell 2 cell 3 cell 4cell 1cell 2cell 3cell 4marker 2marker 1Sitka transformationcell 1cell 2cell 3cell 4Latent marker matrix when ✽✽✽✽ ✽✽ ✽Inferred tree whencell 1 cell 2 cell 3 cell 4latentmarker 2latent marker 1a b c dMarkerPresentAbsentMarkerPresentAbsent1.3. Workflowx = (xc,l), y is the data or a data summary, and ✓ are model parameters. To model errors in copynumber calls, we introduce false positive and negative rate parameters rFP 2 (0, 1) and rFN 2 (0, 1)respectively, and an error matrixerFP,rFN =2641 rFP rFPrFN 1 rFN375 ,p⇣yc,l|xc,l, rFP, rFN⌘= erFP,rFNxc,l,yc,l.from which we set:p(y|x, ✓) =Yl2LYc2Cp⇣yc,l|xc,l, rFPc,l (✓), rFNc,l (✓)⌘.where the scalar erFP,rFNxc,l,yc,l takes values as in Table 1.1. Here we use a global parameterisation wherethe false positive and false negative functions rFPc,l (✓) and rFNc,l (✓) are chosen to be constant.xc,l yc,l pc,l1 0 0 1 rFP2 0 1 rFP3 1 0 rFN4 1 1 1 - rFNTable 1.1: The false positive and false negative errors in sitka’s observation model.In the global parameterisation, we have ✓ = (rFPglobal, rFNglobal). Using a uniform prior distributionfor both error rates can lead to pathological cases as shown in Figure 1.6. To avoid that, we usethe following:rFPglobal ⇠ Uniform⇣0, rFP⌘,rFNglobal ⇠ Uniform⇣0, rFN⌘,We use rFP = 1/10 and rFN = 1/2 as defaults in our experiments. We use a larger rFN; althougherrors in HMM copy number calls can lead to both types of errors, there are more mechanisms thatcan cause a false negative, including a copy number gain, followed by an overlapping copy numberloss.111.3. Workflowx = (xc,l), y is the data or a data summary, and ✓ are model parameters. To model errors in copynumber calls, we introduce false positive and negative rate parameters rFP 2 (0, 1) and rFN 2 (0, 1)respectively, and an error matrixerFP,rFN =2641 rFP rFPrFN 1 rFN375 ,p⇣yc,l|xc,l, rFP, rFN⌘= erFP,rFNxc,l,yc,l.from which we set:p(y|x, ✓) =Yl2LYc2Cp⇣yc,l|xc,l, rFPc,l (✓), rFNc,l (✓)⌘.where the scalar erFP,rFNxc,l,yc,l takes values as in Table 1.1. Here we use a global parameterisation wherethe false positive and false negative functions rFPc,l (✓) and rFNc,l (✓) are chosen to be constant.xc,l yc,l pc,l1 0 0 1 rFP2 0 1 rFP3 1 0 rFN4 1 1 1 - rFNTable 1.1: The false positive and false negative errors in sitka’s observation model.In the global parameterisation, we have ✓ = (rFPglobal, rFNglobal). Usi g a uniform prior distributionfor both error rates can lead to pathological cases as shown in Figure 1.6. To avoid that, we usethe following:rFPglobal ⇠ Uniform⇣0, rFP⌘,rFNglobal ⇠ Uniform⇣0, rFN⌘,We use rFP = 1/10 and rFN = 1/2 as defaults in our experiments. We use a larger rFN; althougherrors in HMM copy number calls can lead to both types of errors, there are more mechanisms thatcan cause a false negative, including a copy number gain, followed by an overlapping copy numberloss.111.3. Workflowx = (xc,l), y is the data or a data summary, and ✓ are model parameters. To model errors in copynumber calls, we introduce false positive and negative rate parameters rFP 2 (0, 1) and rFN 2 (0, 1)respectively, and an error matrixerFP,rFN =2641 rFP rFPrFN 1 rFN375 ,p⇣yc,l|xc,l, rFP, rFN⌘= erFP, FNxc,l,yc,l.from which we set:p(y|x, ✓) =Yl2LYc2Cp⇣yc,l|xc,l, rFPc,l (✓), rFNc,l (✓)⌘.where the scalar erFP,rFNxc,l,yc,l takes values as in Table 1.1. Here we use a global parameterisation wherethe false positive and false negative functions rFPc,l (✓) and rFNc,l (✓) are chosen to be constant.xc,l yc,l pc,l1 0 0 1 rFP2 0 1 rFP3 1 0 rFN4 1 1 1 - rFNTable 1.1: Th lse positive and false negativ errors in sitka’s observation model.In the global parameterisation, we have ✓ = (rFPgl bal, rFNglobal). Using a uniform prior distributionfor both error rates can lead to pathological cases as shown in Figure 1.6. To avoid that, we usethe following:rFPglobal ⇠ Uniform⇣0, rFP⌘,rFNglobal ⇠ Uniform⇣0, rFN⌘,We use rFP = 1/10 and rFN = 1/2 as defaults in our experiments. We use a larger rFN; althougherrors in HMM copy number calls can lead to both types of errors, there are more mechanisms thatcan cause a false negative, including a copy number gain, followed by an overlapping copy numberloss.111.3. Workflowx = (xc,l), y is the data or a data summary, and ✓ are model parameters. To model errors in copynumber calls, we introduce false positive and negative rate parameters rFP 2 (0, 1) and rFN 2 (0, 1)respectively, and an error matrixerFP,rFN =2641 rFP rFPrFN 1 rFN375 ,p⇣yc,l|xc,l, rFP, rFN⌘= erFP,rFNxc,l,yc,l.from which we set:p(y|x, ✓) =Yl2LYc2Cp⇣yc,l|xc,l, rFPc,l (✓), rFNc,l (✓)⌘.where the scalar erFP,rFNxc,l,yc,l takes values as in Table 1.1. Here we use a global parameterisation wherethe false positive and false negative functions rFPc,l (✓) and rFNc,l (✓) are chosen to be c stant.xc,l yc,l pc,l1 0 0 1 rFP2 0 1 rFP3 1 0 rFN4 1 1 1 - rFNTable 1.1: The false positive and false negative errors in sitka’s observation model.In the global parameterisation, we have ✓ = (rFPglobal, rFNglobal). Using a uniform prior distributionfor both error rates can lead to pathological cases as shown in Figure 1.6. To avoid that, we usethe following:rFPglobal ⇠ Uniform⇣0, rFP⌘,rFNglobal ⇠ Uniform⇣0, rFN⌘,We use rFP = 1/10 and rFN = 1/2 as efaults in our experime ts. We use a larger rFN; althougherrors in HMM copy number calls can lead to both types of errors, there are more mechanis s thatcan cause a false egative, including a copy n mber gain, followed by an overlapping copy numberloss.11,1.3. Workflowx = (xc,l), is the data or a data summary, and ✓ ar model parameters. To model errors in copynumber c lls, we introduc false positive and negative rate parameters rFP 2 (0, 1) and rFN 2 (0, 1)respectively, and an error matrixerFP,rFN =2641 rFP rFPrFN 1 rFN375 ,p⇣yc,l|xc,l, rFP, rFN⌘= erFP,rFNxc,l,yc,l.from which we set:p y|x, ✓) =Y2LYc2Cp⇣yc,l|xc,l, rFPc,l (✓), rFNc,l (✓)⌘.where the scala erFP,rFNxc,l,yc,l takes values as in T ble 1.1. Here we use a global parameterisation wherethe lse p sitive and fals negative functio s rFPc,l (✓) and rFNc,l (✓) are chosen to be constant.xc,l yc,l pc,l1 0 0 1 rFP2 0 1 rFP3 1 0 rFN4 1 1 1 - rFNTable 1.1: The false positive and fal negative erro s in sitka’s observation model.In t e global parameterisation, we have ✓ = (rFPglobal, rFNglobal). Using a uniform prior distributionfor both error rates can lead to pathological cases as s own in Figure 1.6. To avoid that, we usethe following:rFPglobal ⇠ Uniform⇣0, rFP⌘,rFNglobal ⇠ Uniform⇣0, rFN⌘,We use rFP = 1/10 and rFN = 1/2 as defaults in our experiments. We use a larger rFN; althougherrors in HMM copy number call can l ad to both types of errors, there are more mechanisms thatca c use a false n gative, including a copy number ain, followed by an overlapping copy numberloss.111.3. Workflowx = (xc,l), y is the data or a data sum ary, and ✓ are model parameters. To model errors in copynumber calls, we introduce false positive and negative rate parameters rFP 2 (0, 1) and rFN 2 (0, 1)respectively, and an error matrixerFP,rFN =2641 rFP rFPrFN 1 rFN375 ,p⇣yc,l|xc,l, rFP, rFN⌘= erFP,rFNxc,l,yc,l.from which we set:p(y|x, ✓) =Yl2LYc2Cp⇣yc,l|xc,l, rFPc,l (✓), rFNc,l (✓)⌘.where the scalar erFP,rFNxc,l,yc,l t kes v lues as in Table 1.1. H we use a global parameterisation wherethe false positive and false negative functions rFPc,l (✓) and rFNc,l (✓) are chosen to be constant.xc,l yc,l pc,l1 0 0 1 rFP2 0 1 rFP3 1 0 rFN4 1 1 1 - rFNTable 1.1: The false positive and false negative errors in sitka’s observation model.In the global parameterisation, we have ✓ = (rFPglobal, rFNglobal). Usi g a uniform prior distributionfor both error rates can lead to pathological cases as shown in Figure 1.6. To avoid that, we usethe following:rFPglobal ⇠ Uniform⇣0, rFP⌘,rFNglobal ⇠ Uniform⇣0, rFN⌘,We use rFP = 1/10 a d rFN = 1/2 as d faults in our experiments. We use a larger rFN; althougherrors in HMM c py number calls can lead to both types of errors, there are more mechanisms thatcan cause a false negative, including a copy number gain, followed by an overlapping copy numberloss.111.3. Workflowx = (xc,l), y is the data or a data summary, an ✓ a e model arameters. To model errors in copynumber calls, we int oduce f ls po itive and neg tive rate parameters rFP 2 (0, 1) and rFN 2 (0, 1)respectively, and n error matrixerFP,rFN =2641 rFP rFPrFN 1 rFN375 ,p⇣yc,l|xc,l, rFP, rFN⌘= erFP, FNxc,l,yc,l.from which we set:p(y|x, ✓) =Yl2LYc2Cp⇣yc,l|xc,l, rFPc,l (✓), rFNc,l (✓)⌘.where the scalar erFP,rFNxc,l,yc,l take v lues as in Table 1.1. Here we use a global parameterisation wherethe false positive and false negative functions rFPc,l (✓) d rFNc,l (✓) are chosen to be constant.xc,l ,l pc,l1 0 1 rFP2 0 1 rFP3 1 0 rFN4 1 1 1 - rFNTable 1.1: Th false po tive and false negative rrors in sitka’s observation model.In the global parameterisation, we have ✓ = (rFPgl bal,FNglobal). Using a uniform prior distributionfor b th rror rate can lead t pathological cases as sho n in Figure 1.6. To avoid that, we usethe following:rFPglobal ⇠ Uniform⇣0, rFP⌘,rFNglobal ⇠ Uniform⇣0, rFN⌘,We use rFP = 1/10 and rFN = 1/2 as defaults in our exp ime ts. We se a larger rFN; althougherrors in HMM copy number calls can lead t both types of errors, there are more mechanisms thatcan cause a fals negative, including a copy number gain, followed by an overlapping copy numberloss.111.3. Workflowx = (xc,l), y is the data or a data sum ry, and ✓ are mod l parameters. To model errors in copynumb r calls, we introduce false positive and negative rate parameters rFP 2 (0, 1) and rFN 2 (0, 1)respectively, and an error matrixerFP,rFN =2641 rFP rFPrFN 1 rFN375 ,p⇣c,l|xc,l, rFP, rFN⌘= erFP,rFNxc,l,yc,l.from which we set:p(y|x, ✓) =Yl2LYc2Cp⇣yc,l|xc,l, rFPc,l (✓), rFNc,l (✓)⌘.where the scalar erFP,rFNxc,l,yc,l t kes v lues a in Ta le 1.1. H we use a global parameterisation wherethe false positive and false negative functions rFPc,l (✓) and rFNc,l (✓) are chosen to be c stant.xc,l yc,l pc,l1 0 0 1 rFP2 0 1 rFP3 1 0 rFN4 1 1 1 - rFNTable 1.1: The false positive a d false negative errors in sitka’s observation model.In the glob l parameterisation, we have ✓ = (rFPglobal, rFNglobal). Using a uniform prior distributionfor both err r rates n le d to pathological cases as shown in Figure 1.6. To avoid that, we usethe following:rFPglobal ⇠ Uniform⇣0, rFP⌘,rFNglobal ⇠ Uniform⇣0, rFN⌘,We use rFP = 1/10 and rFN = 1/2 as faults in our experime ts. We use a larger rFN; althougherrors in HMM c py number calls can lead to both types of errors, there are more mechanis s thatca cause a false egativ , including a copy n mber gain, foll wed by an overlapping copy numberloss.11,marker 2marker 1root rootFigure 2.7: Pathological tree reconstruction under default observation prior. (a) The true treereconstruction in a simple example with a balanced phylogeny with two clades of size two, andtwo unique markers, coloured red and blue, that distinguish the left and right clades respectively.(b) The binarised input matrix corresponding to the four cells at the two markers. The desiredobservation error rates should be zero and the latent and observed marker matrices should matchexactly, as the perfect phylogeny assumption holds. If the observation error parameters are set toone, that is rFPglobal = 1 and rFNglobal = 1, then the latent marker matrix with all entries flipped asshown in (c) will have an equal likelihood under this setting as the desired latent matrix has whenerror rates are set to zero (see Equation (2.4)). (d) The incorrect tree reconstruction where the leftand right clades are erroneously assigned to the blue and red markers.We represent the prior on perfect phylogenetic trees via the following two step generatingprocess: (i) sample a mutation tree topology t (i.e., a type I tree with all cell vertices removed),(ii) assign all cells c to t. Recall that the tree t comprises all markers l ∈ L as well as an additionalimaginary marker v∗ that is the root. If a marker l′ is an ancestor of marker l in t, we assume thatall cells that are descendants of l have the trait l′ as well. If a cell c is attached to marker l, weposit that it has a copy number change at all the markers in the shortest path from l to the rootmarker v∗ and has no copy number changes across the other markers in L.For simplicity, we use a uniform prior on tree topologies as follows:p(t) =1 [t ∈ T ](|L|+ 1)|L|+|C|−1 ,where T is the set of all perfect phylogenetic trees that result from the two step generative process24described above. Simulation from the prior can be performed using Wilson’s algorithm [57], followedby independent categorical sampling to simulate the cell assignments.2.3.5 InferenceWe aim to approximate the posterior distribution of the tree t and the parameters θ via MCMCpi(t, θ) ∝ p(t)p(θ)p(y|x(t), θ).First we describe a tree exploration MCMC move that considers an exponential number of neigh-bouring trees, then we review how parameters rFPglobal and rFNglobal are sampled.sitka uses a tree sampling move to explore a large neighbourhood of a given tree. Given a treet and marker l, we remove l from t, and add a new marker such that the new tree t′ remains aperfect phylogeny, i.e., t′ ∈ N l(t) ⊂ T . After removing l (edge contraction), adding a new markerthat is an edge insertion move can be described as follows:1. Pick a non-cell vertex v from R = {v∗} ∪ L\{l} where v∗ is the root node.2. Pick any subset of v’s descendent subtrees and remove it from v.3. Add a new node v′ under v and move the selected nodes from step 2 above and attach themto v′.Note that the choice of node v in the first step above, partitions N l(t) into blocks of N lv(t\l),that is N l(t\l) = ∪vN lv(t\l). Figure 2.4-e (right) shows an example of an edge insertion move. Amarker named chr15 5950, coloured red, has three children, all cell vertices in the sampled tree atMCMC iteration 100. This would be node v. In the next iteration, two of its children, namelycells RC07C and RC05C4 are chosen and removed from v. They are then inserted under markerchr1 4900, v′, which is now a child of marker chr15 5950.The probabilities required in step 1 above are of the form:ρ¯v =ρv∑v˜∈R ρv˜,25where:ρv =∑t∈N lv(t\l)p(t)p(y |x(t), θ), (2.5)and t\l denotes the tree t with node l removed. To compute ρv, we start with the following recursionfor all vertices v in t\l: first, for all vertices c corresponding to a cell and b ∈ {0, 1}, define:pbc = p(yc,l | b, θ),where p(yc,l | b, θ)is defined in Equation (2.4). Next, we compute these distributions for all subtreesof t\l. This can be done efficiently via a bottom-up recursion on the tree t\l: for all v ∈ R, b ∈ {0, 1},pbv =∏v′′∈ children(v)pbv′′ ,where children(v) denotes the list of children of vertex v.It can be shown that we can compute the probabilities required in step 1 above as follows ([1]):ρ¯v =ρv∑v˜∈R ρv˜(2.6)=(∏vi∈children(v)(p0vi+p1vi)p0v)∑v˜∈R∏v′i∈children(v˜)(p0v′i+p1v′i)p0v˜. (2.7)Once v is sampled, we choose a subset of its children to move to v′ by sampling k independentBernoulli random variables with the probability for each child vi of v as below and selecting childrenwith corresponding Bernoulli realisations of 1:p1vip0vi + p1vi.To resample the parameters θ we use a slice sampling algorithm [58] and condition on the26hidden state matrix x in a Metropolis-within-Gibbs framework. Briefly, we compute two sufficientstatistics from the matrix x, (i) the number of false positive instances, nFP, and (ii) the number offalse negative instances, nFN as follows:nFP = nFP(x) =∑c∈C∑l∈L1[xc,l = 0, yc,l = 1]nFN = nFN(x) =∑c∈C∑l∈L1[xc,l = 1, yc,l = 0].Based on these sufficient statistics, we obtain:p(y |x, θglobal) =(rFPglobal)nFP(rFNglobal)nFN(1− rFPglobal)nA−nFN(1− rFNglobal)nP−nFP, (2.8)where nP and nA denote the number of present and absent markers in the observed data and canbe pre-computed as follows:nP =∑c∈C∑l∈L1[yc,l = 1]nA = |C||L| − nP.The model is implemented in the Blang probabilistic programming language [59]. We use aparallel tempering (PT) algorithm for inference (see [1] for a detailed explanation of the annealeddistributions) and initialise it via a tree sampled uniformly from the prior.2.3.6 Summarising the posterior distributionHere we approximate the Bayes estimator by minimising the Bayes risk [60]:argmint∈T∑t′∈T∫L(t, t′)pi(t, dθ). (2.9)27where pi(t, θ) is the posterior distribution of the tree t and the parameters θ and L(t, t′) is a lossfunction. We use the L1 metric on the matrices of induced indicators x(t) as the loss-function:L(t, t′) =∑l∈L∑c∈C|xc,l(t)− xc,l(t′)|.It is useful to define marginals mc,l that could be conceptualised as the posterior probability ofcell c to have trait l:mc,l =∑t∈T∫1[xc,l(t) = 1]pi(t, dθ),Using the MCMC samples t1, t2, . . . , tN to average marginals mc,l we get a Monte Carlo ap-proximation:1NN∑i=1xc,l(ti)→ mc,l a.s.Figure 2.4-g shows an example of the matrixm, each element of which is one of the approximatedm¯c,l . We can now write the objective function of Equation (2.9) via the above marginals:∑t′∈T∫L(t, t′)pi(t, dθ) =∑t′∈T∫ ∑l∈L∑c∈C|xc,l(t)− xc,l(t′)|pi(t, dθ)=∑l∈L∑c∈C∑t′∈T∫|xc,l(t)− xc,l(t′)|pi(t, dθ)=∑l∈L∑c∈C{mc,l(1− xc,l(t)) + (1−mc,l)xc,l(t)}=∑l∈L∑c∈C{xc,l(t)− 2mc,lxc,l(t)}+ constant (2.10)We use a greedy algorithm to approximately minimise the Equation (2.10). Briefly, we startwith a star-shaped tree with leaves C rooted at v∗ and add markers from L one by one from amarker queue sorted by priority score. The priority score of each marker l is computed aspriority(l) = maxt′∈N l(t)q(t′)∑t′′∈N l(t) q(t′′)where28q(x) =∏c∈C∏l∈L(x)qc,l(xc,l)qc,l(xc,l) = 2mc,lxc,l − xc,l.The quantities in the priority queue can be computed as in Section 2.3.5. We take the result ofminimising the Bayes risk as the consensus tree.2.3.7 Consensus tree and CNA heatmap visualisationIn visualising the output of sitka, we always use the type II trees. In what follows we introducetwo different ways of visualising sitka trees that associate the tree to their input and parameterspace.Marker-centric viewThis view maps the location of a marker to its position on the genome. As each marker is associatedwith a pair of consecutive bins, we take the convention of setting the genomic position of a marker tothat of its associated leftmost bin. The plots comprise a type II tree t (Section 2.2) on the left nextto a (L−1)×(L−1) matrix w to the right. In this plot, each node on the tree corresponds to a row.The matrix w has exactly L−1 non-empty elements. Each element wi,j is non-empty if marker i ison the i-th row of t and the j-th bin. The direction of the change in CN corresponding to a markeris shown by a left-pointing or right-pointing arrow for increase or decrease in CN respectively.Figure 2.8 shows the marker-centric view for the SA535 dataset. The rows of the matrix can beannotated to show information about specific internal nodes. For instance, in Figure 2.8, the Alland Immed columns show what fraction of the genome is altered between a node and its parenton the tree. The former constructs a summary genotype for the parent node based on cells thatare directly attached to it (if they exist), while the latter uses all cells that are descendent of theparent node to do so. We use the CNA profiles to compute the summary genotype. Let asub be aninteger valued Csub ×Dsub matrix that encodes the CN state of a set of Csub cells over Dsub copynumber bins. We define s as a column vector of size Dsub where si = fsummary(a1:Csub,i). That29is, s is an artificial cell whose CNA profile reflects that of those in asub. We set fsummary(.) to themedian function. To compute the percentage of the genome altered, we compute the Manhattandistance between the median genotypes of each node and its parent, normalised by the number ofcopy number bins. We acknowledge that the summary genotype at an internal node computed asabove is a heuristic for data exploration purposes only as the phylogenetic tree is based on the CNchanges, and not the CNs. In contract, the marker presence or absence profile can be reconstructedat any internal node v via the function x(v).30| | | | | | | | | | | | | | | | | | | | || | | | | | | | | | | | | | | | | | |1 2 3 4 5 6 7 8 9 10 11 12 13 15 16 17 18 19 20 21 X►◄◄►◄◄◄◄◄◄◄◄►◄◄►◄►◄►◄►◄►►◄►►◄◄◄◄►►►►►◄◄◄►◄►◄►◄►◄►◄◄►►►►►◄►►◄►◄►►►►►◄►◄►►◄◄►◄►◄◄►►►►►◄►►►►►►►►◄►►◄◄◄►►◄◄◄◄◄►►◄◄◄◄◄◄►◄►►►►◄►►◄►◄►►►►►►◄►►►| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | |AllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllAllImmed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Immed.Cell counts 5 10 15 200.5Averages0.035 0.040 0.045 0.050 0.055Fraction of genome altered versus parent (per bin) Breakpoint sign  Positive Negative#Cells = 493 -- #Markers = 139SA535AllImmediateMarker Figure 2.8: Marker-centric view for SA535. A type II tree where all cell nodes are dropped andchains are collapsed. The terminal nodes on the tree represent markers and their colours map togenomic positions as indicated on top of the chromosome (horizontal) axis (e.g., yellow maps tochromosome 2). Note that there are more negative (right-facing triangle) than positive (left-facingtriangle) CN changes on chromosome 12 (7 and 4 respectively). One possible explanation is that oneend of at least 3 CN changes are unobserved; potentially the imbalance may be due to sequencingerrors. Another explanation is that the change points are located or stretch over to the end of thechromosome and therefore are not reflected in sitka’s lossy binary transformation.31Multi-channel viewOne way to visualise the tree inference results is to arrange the tree and the cell-by-locus CNmatrix side by side where the rows of the matrix correspond to the position of individual cellson the tree and the loci are arranged by their genomic position. This is shown in Figure 2.4-h.Figure 2.12-a-c are examples of the multi-channel visualisation where each marker is representedby a tuple of three different data-types, namely, (i) the markers, (ii) the sitka posterior, and(iii) the sitka transformation. We describe each channel for a fixed marker l in tree t. The firstchannel corresponds to the induced matrix x(t) and indicates which cells in the consensus tree aredescendants of l. In other words, whether or not in the latent matrix, the inferred marker is presentin each cell. The second channel corresponds to the marginal matrix m¯ as defined in Section 2.3.6,and shown in Figure 2.4-g. It denotes the posterior probability of the trait l to be present in eachcell. Finally, the third channel corresponds to the transformed marker matrix y (see Section 2.3.3)and shows whether the trait was observed in each cell.This view can be used to quickly assess the discrepancy between the input data and the inferredtree. It may illuminate cases in which the perfect phylogeny assumption is violated. An exampleis ChrX in the OV 2295 dataset (Figure 2.12-a). ChrX has a long orange band (inferred marker)not matched by a black band (jitter-fixed observed marker) suggesting that a perfect phylogenyviolation may have occurred. The pattern in this marker is consistent with the presence of anancestral event followed by a deletion.In Figure 2.12-b, a set of diploid cells are attached to the root of the tree. These are control cellsincluded in the experiment and correspond to a marker region in the bottom of the matrix with noinferred markers (orange bands) and almost no observed markers (black bands). In this dataset,there are multiple change points where the observed marker has a high density (black band), butthe tree is reconstructed with the marker absent (no matching orange band). Examples includeChr1, Chr7 and Chr16. One possible explanation could be that the end-points of each event weredetected as slightly shifted in each cell. For instance, in Figure 2.11 there are two loci with anamplification (CN state equal to three) in the p-arm of Chr1 where the first locus is amplified infewer cells than the second locus. In particular, cells that harbour a mutation in the first locus,appear to not have a mutation in the second locus, suggesting that the same event was called in32the first locus in some cells, and in the second locus in others. A solution may be to fix jitter with alarger neighbourhood size (Section 2.3.3). An alternative hypothesis is that the cells in this datasethave a mutator phenotype that promotes de novo CN mutations in these loci.Another locus suspected of violating the perfect phylogeny assumption is in the SA535 dataseton Chr3 (Figure 2.12-c). A black band is not matched by an orange one. This pattern could bea result of two insertion events, or two deletion events. For instance, in the former, the markeris hypothesised to be absent in the ancestral state, and it is gained in two subsequent clades (thesubset of the tree aligned with the black band). While in the latter, the converse holds.Define the mismatch rate as the fraction of cells that have a discrepancy between the jitter fixedinput and the latent tree. Figure A.4 shows the distribution of mismatch rate for each dataset.The distribution for OV 2295 and SA535 have an approximately bimodal shape while the SA501distribution exhibits a heavy tail. In OV 2295, 41 markers (11%) have a mismatch rate of over50%, where marker chr15 67000001 67500000 has the highest mismatch rate at 70%. In SA501,30 markers (11%) have a mismatch rate of over 50%, 13 of which (5%) have a mismatch rate ofover 75%. SA535 has the lowest maximum mismatch rate at 49% (marker 15 72000001 72500000 ).While in the OV 2295 dataset, and to a lesser extent, the SA535 dataset, there is some evidence of acluster of markers with a high mismatch rate, the distribution is not concentrated at the minimumand maximum mismatch rates. This suggests that a few markers do not explain a large portion ofthe observed noise.2.4 Overview of datasetsHere we introduce three real-world datasets and look at their type II trees. The first dataset,SA535 from [20], contains 679 cells from three passages of a triple negative breast cancer (TNBC)patient derived xenograft sample. Passages X1, X5, and X8 had 62, 369, and 231 cells post qualityfiltering respectively. We also included 17 mostly diploid control cells. These cells were combined togenerate the input to the analysis pipeline. Figure 2.9 shows the phylogenetic tree and the heatmapencoding of the CNA profiles where the cells are sorted by their order appearance on the tree. Thesecond dataset labelled OV A consists of cells from three samples taken from a patient with highgrade serous (HGS) ovarian cancer; the first sample, SA1090, was from an ascites pre-treatment,33while SA922 was from an ascites post-treatment. The third sample, SA921, was taken from theovary. See Figure 2.10 for the tree and the CNA profile heatmap for this dataset. The final dataset,SA501 [43], is another TNBC xenograft tumour from 6 untreated passages, namely X2, X5, X6,X8, X11, and X15. After filtering we had 515, 236, 328, 189, 836, and 308 cells in each passagerespectively (for a total of 2,412 cells, see Figure 2.11). Table A.5 shows the attrition after eachstep of filtering cells per passage in each dataset. We discuss how CNA rate changes over time withthe passage information in Section A.5.Figure 2.9: Phylogenetic tree and CNA profile heatmap for the SA535 dataset. The rows of theheatmap are sorted according to the placement of cells on the phylogenetic tree. The columns ofthe heatmap are sorted by their genomic position.34Figure 2.10: Phylogenetic tree and CNA profile heatmap for the OV A dataset. The nearly diploidcells with the loss of heterozygosity on chromosome X are from SA1090. The cells with an ampli-fication on chromosome 22 are from SA922. The rest belong to SA921.Figure 2.11: Phylogenetic tree and CNA profile heatmap for the SA501 dataset. Note that thediploid cells at the bottom of the heatmap are control cells that were included in the experiment.2.5 Tree evaluation2.5.1 Predictive testTo evaluate the inferred trees, we suggest a test that involves predicting the entries in the inputbinary matrix given to the tree inference method. We take the binarised input matrix y, the inputmatrix to the sitka algorithm as described in Section 2.3.3 as ground truth. Consider an inferred35tree, t, and the corresponding matrix g = x(t). In general the inferred trees from the baselinemethods do not have named internal nodes, nor do they have the same number of internal nodes asthe number of marker L. Therefore we do not know which marker in the inferred tree t correspondsto which column in the matrix y. We note that this is not the case with trees inferred from sitkawhere the internal nodes of the tree t correspond to the columns of the induced matrix x(t). As aresult, for methods other than sitka, for each column in the input data matrix, we pick a clade int that has the highest prediction accuracy for the entries in that column.For each method, we report Youden’s J index [61]. We define below the function h to be abinary classification counts matrix, i.e., for two column vectors w and z of size C, it forms theconfusion matrix. h : {0, 1}C × {0, 1}C → {0, 1}2×2 wherehi,j(w, z) =∑c∈C1 [wc = i, zc = j] .For example h0,0(w, z) would count the number of times both elements of w and z were equal tozero (or the number of true negatives). We define accuracy for a given confusion matrix ζ computedfrom the h map above as:acc(ζ) :=ζ0,0 + ζ1,1∑i∈{0,1}∑j∈{0,1} ζi,jWe further define sensitivity and specificity assensitivity(ζ) :=ζ1,1ζ1,1 + ζ1,0specificity(ζ) :=ζ0,0ζ0,0 + ζ0,1youden(ζ) := sensitivity(ζ) + specificity(ζ)− 1For a given tree t and its latent matrix g = x(t) we compute the Youden’s index as follows:1. For all marker l in y, compute ζl = argmaxζl′ ,l′∈columns(g) acc(ζl′).2. For tree t, set compute the confusion matrix as ζt =∑l′∈columns(g) ζl′ .3. Define the Youden’s index for t as youdent := youden(ζt).That is for each marker in y, we take the clade that among all possible clades in t maximises36the accuracy in predicting which cells are present in the l-th column of y. We then sum over allthese scores to compute a confusion matrix for t and use this agglomerative matrix to compute theYouden’s index for the tree. We use the delta method to calculate confidence intervals. Figure 2.12-d shows the Youden’s index, its 95% confidence interval for sitka and 6 baseline methods over 3different real-world datasets. sitka has a higher score than all competing methods. In Section 2.5.2we review these results in detail.A similar idea can be used to compare the different binarisation methods discussed in Sec-tion 2.3.2. Briefly, we report a leave-one-out score that represents how well the model predictsy(c,l) | y\(c,l) where y\(c,l) denotes all the elements of the matrix except one entry. In other words,how well the model predicts one entry of the y matrix given all the other entries. This is averagedover all cell-marker pairs (c, l). It is computed using importance sampling (in the context of efficientleave-one-out estimation [62]). To compare different encodings of the input matrix, a common setof markers should be used. For example to compare preprocessing methods 1 and 2 that resultin sets of markers L1 and L2, respectively, we ensure L1 ∩ L2 6= ∅, then the leave-one-out score isaveraged over this intersection only.2.5.2 Comparison to baseline methodsHere we compare the performance of sitka against a number of baseline methods over real-worlddatasets. We first briefly overview the baseline methods. The comparisons are summarised inFigure 2.12-d.Overview of baseline methodsHere we briefly review the baseline phylogenetic methods. For all agglomerative clustering methods,the distance function d(c1, c2) between two single-cells c1 and c2 was set to the Euclidean distancebetween their CNA profiles.UPGMA (Unweighted Pair-Group Method with Arithmetic mean) is an agglomerative clusteringmethod that, given a similarity matrix, constructs a rooted phylogenetic tree [63]. We start byputting each datapoint in a separate cluster. At each step of the algorithm, we merge the twoclusters with minimum inter-cluster distance. Inter-cluster distance between clusters X and Y is37defined by:d (X,Y ) =1|X||Y |∑x∈X,y∈Yd (x, y) (2.11)where |.| denotes the cardinality of a set. In other words, the average pairwise distance betweenelements of X and Y .WPGMA (Weighted Pair-Group Method with Averaging) is similar to UPGMA but uses the followingto compute the distance between clusters X and Y as defined by [63]:d (X,Y ) =12∑x∈X,y∈Yd (x, y) (2.12)In other words, it ignores the cardinality of clusters when computing the distance between them.NJ (Neighbour joining) is a similar agglomerative clustering method that accepts an initialdistance matrix d [64]. It keeps a running distance matrix q between each node of the tree andupdates it over each iteration. In each iteration, NJ (i) re-calculates q, (ii) picks (f, g) to joinunder a new node u such that f 6= g and q(f, g) is the minimum entry in q, and (iii) calculatesbranch-length for the two newly created edges (u, f) and (u, g). The running distance matrix q isupdated as follows:q(i, j) = (n− 2)d(i, j)−n∑k=1d(i, k)−n∑k=1d (j, k)The entries for nodes (f, g) in matrix d is then replaced by that of the newly created node u wherefor each node k /∈ {f, g} we have:d(u, k) =12(d(f, k) + d(f, g)− d(f, g))HDBSCAN first computes a 2-dimensional representation of the copy number matrix, i.e., caststhe C × L CNA matrix a into a C × 2 embedding matrix u. This step uses UMAP [65]. Secondly,a hierarchical clustering is computed via the method described in [66].MrBayes is a Bayesian phylogenetics framework that implements multiple evolutionary modelsand uses MCMC to approximate the posterior distribution of trees and model parameters [67].38We ran MrBayes version 3.2.6 on two types of copy number input matrices: (i) integer copy num-bers from 0 to 9 across the entire genome where CNs greater than 9 were assigned a value of 9(MrBayes-np2 or MrBayes-np8), and (ii) the same as sitka’s input of binary copy number changesat a reduced set of genomic locations (MrBayesWithBinaryInput-np2 or MrBayesWithBinaryInput-np8).We ran MrBayes using either two runs with one chain each (MrBayes-np2, MrBayesWithBinaryInput-np2),or two runs with four chains each (MrBayes-np8 or MrBayesWithBinaryInput-n8).MrBayesWithBinaryInput-np8 did not finish running on any of our datasets. This is expectedas the model based assessment of likelihood usually leads to undesirably long execution time. Insome cases it is possible to speed up the run time of MrBayes via the use of graphical process-ing units (GPUs) and parallelization ([68]). However, it is unlikely to scale up to the sizes ofdatasets that we have in this work, and anticipate would be soon available. One reason is thatthe parallelization occurs across sites (markers) and not taxa (cells). For brevity, we will denoteMrBayesWithBinaryInput-np2 by MrBayesWithBinaryInput. Following [43], we ran 10 millioniterations with 50 percent burn-in fraction. All parameters of the likelihood model were left at thedefault values except lset rates=adgamma in order to account for correlated rates for adjacentsites. When using integer CN states as input to MrBayes, we set datatype=standard.Benchmarking results using the predictive testNow we use the predictive test derived in Section 2.5.1 to compare sitka to baseline methodsreviewed above, over three real-world datasets SA535, OV A, and SA535. Figure 2.12-d shows theYouden’s J index and its 95% confidence interval for each method.sitka has the highest Youden’s index across all three datasets. UPGMA and WPGMA perform sim-ilarly on SA501 and SA535. UPGMA does slightly better than WPGMA on the OV A dataset. HDBSCANhas a close but slightly smaller Youden’s index than UPGMA over the SA535 and OV A datasets,but performs marginally better on SA501. NJ trails WPGMA on SA501 and the OV A datasets, andhas the lowest Youden’s index on SA535. MrBayes does well on the smallest dataset, SA535, withMrBayes-np2 and MrBayes-np8 performing similar to WPGMA, and MrBayesWithBinaryInput hav-ing achieved the second highest Youden’s index. On the OV A data, MrBayesWithBinaryInputand MrBayes-np2 trail behind NJ, while MrBayes-np2 has the lowest Youden’s index among allmethods on all datasets. Similar to the OV A case, MrBayesWithBinaryInput and MrBayes-np239trail behind NJ over the SA501 dataset. MrBayes-np8 did not finish running on SA501 after severaldays. This suggests that the algorithm used in MrBayes needs prohibitively more computationalbudget to achieve a Youden’s index on par with the other methods. The results in this comparisonsuggest that sitka performs better than the baseline methods and therefore we chose it to inferclonal phylogenies in Chapter 4.tree_mat ● ●0 10.00 0.25 0.50 0.75 1.00average filtered● ●0 1| | | | | | | | | | | | | | | | | | | | | | || | | | | | | | | | | | | | | | | | | | | | |1 2 3 4 5 6 7 8 9 10111213141516171819202122 X| ||| ||| | | || ||| ||| || |||| ||| | | | || || | |||| || || | ||||| |||| || | | || | || | || || | || || | || ||| || ||| || || || |||||| |||| ||| || || | | || | | | |||| ||| | || ||| || | |||| || || | || | ||||| | ||| | | || | | |OVAOVA, 997 cells, 360 markerstree_mat ● ●0 10.00 0.25 0.50 0.75 1.00average filtered● ●0 1| | | | | | | | | | | | | | | | | | | | | || | | | | | | | | | | | | | | | | | | | | |1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22|| | || || | || || | || || | ||| || ||||| || || | | ||| | | | | || | ||| | | | || ||| | | || | ||| || || | ||| || | | |||| || || ||||| | || || | | || ||| ||| || || | || | || || || ||| | | || || | | || | ||| |||| | | ||| ||| ||| | ||SA501CellsSA501, 2412 cells, 263 markersCellsMarkers Markerstree_mat ● ●0 10.00 0.25 0.50 0.75 1.00average filtered● ●0 1| | | | | | | | | | | | | | | | | | | | | | || | | | | | | | | | | | | | | | | | | | | | |1 2 3 4 5 6 7 8 9 10111213141516171819202122 X| || || ||| || | |||| |||| || |||| |||| ||| || ||| || | || | || || || ||| ||| |||||| || | || | || ||| | || | ||| | | || | | || | |||| || ||| || || | ||| ||| | ||| | | || || | || || || || |SA535CellsMarkersadcbSA535, 493 cells, 391 markers Sitka transformationPresentAbsentMarkersPresentAbsentSitka posteriorPresent0.0 1.0Benchmarking versus baseline methods●●●●●●●●●●●● ●●●●●●● ●●● ●0.50.60.70.80.91.0SA501 TNBC−SA535 OVA−2295DatasetYouden's indexTree inference method ●●●●●●●●sitkaUPGMAHDBSCANWPGMANJMrBayes−np2MrBayes−np8MrBayesWithBinaryInput−np2Sitka transformationPresentAbsentMarkersPresentAbsentSitka posteriorPresent0.0 1.0Sitka transformationPresentAbsentMarkersPresentAbsentSitka posteriorPresent0.0 1.0Youden’s index95% confidence intervalFigure 2.12: Results over real-world datasets and benchmarking against baseline methods. (a), (b),and (c) show the consensus tree and the multi-channel marker matrix for the OVA, SA501, andSA535 datasets respectively. In each panel, the tree is the consensus tree constructed as describedin Section 2.3.6 and the matrix is the multi-channel visualisation of Section 2.3.7. Each markeris represented by a marker tuple, comprising three consecutive columns we call channels (as inFigure 2.4-i). The three channels encode for each marker in each cell (i) whether it is inferred tobe present, (ii) the posterior probability of its presence, and (iii) whether it was present in thetransformed input matrix. The number of rows (cells) and markers (columns) are noted on top ofeach panel. The position of the internal nodes on the genome is colour-coded to match the colouron the chromosome axis. (d) Comparison of methods. The dot shows the Youden’s J index, whilethe bar represents the 95% confidence interval.402.6 Identifying clones through phylogenetic analysisUsing sitka we can establish the evolutionary relationships of cells in a heterogeneous sample. Toinvestigate cancer evolution we need to determine the abundance of subpopulations over time. Tothis end, we introduce Lumberjack, a tree-cutting algorithm that we use to define clonal subpop-ulations. In the output tree of sitka, cells are part of the terminal leaf nodes of the phylogenetictopology. We post-process the inferred trees to identify clonal populations from major clades.When clonal populations are defined, their abundances can be counted as a function of timeseriesand these can be used for fitness inference (see Chapter 3). Clones are constructed by identifyingconnected components (each a clade or a paraphyly) in the phylogenetic tree reconstruction. Thetree is ‘cut’ into discrete populations according to the following procedure (see Algorithm (2)).We review and build on the notation used in this chapter to facilitate the introduction of theLumberjack algorithm. Let L be a set of markers and C be a set of cells. Define t = (L,C,E) to bea rooted phylogenetic tree with E its set of directed edges. We assume that t is a type II tree (Sec-tion 2.2). Let |t| = |C|, that is the number of cells that belong to tree t. Let tl = (Ll, Cl, El) denotethe subtree rooted at node l. Define parent(l) to be the parent of node l, and let descendants(l)comprise all its descendants.The inputs to the algorithm are the rooted phylogenetic tree t, the CN states of its cells andthe minimum Mmin and maximum Mmax allowed clone sizes. A clone is defined as a connectedcomponent (each a clade or a paraphyly) in the graph tree t composed of cells of sufficient genomichomogeneity. The degree of homogeneity can be tuned by limiting the number of markers and thedifference in the CN of sub-clades in a clone. The algorithm first finds the coarse structure, thatis, it divides the tree into major clades, and then looks for fine structures within each clade bytraversing the tree in a bottom up manner and merging markers that are sufficiently similar. Theremaining markers constitute the roots of detected clades (clones).To obtain the coarse structure from the reconstructed phylogenetic tree we use a two stepprocedure: (i) identify monophyletic clades via Algorithm (2), (ii) then remove the cells comprisingthe clades found in step one from the tree and repeat Algorithm (2). We note that these new clades(if any) could be paraphyletic.To find the fine structures within the initial clades we use the following procedure. For each clade41Algorithm 2 Heuristic (Top-down)1: procedure sitka-lumberjack(t, Mmin, Mmax)2: marker-queue ← depthFirstSearch(t)3: for marker l in marker-queue do4: if m ≤ |tl| ≤M then5: CUTS← CUTS ∪ l6: Remove all markers below l from t7: else8: for marker l′ in descendants(l) do9: if Mmin ≤ |tl′ | ≤Mmax then10: CUTS← CUTS ∪ l′11: Remove all markers below l′ from t12: if no eligible markers were found then13: CUTS← CUTS ∪ l14: Remove tl from t15: return CUTSρ, and its corresponding sub-tree tρ, denote by Lρ a set of markers l for which Mmin ≤ |tl| ≤Mmax.In a bottom-up traverse of the tree, for each node l ∈ Lρ, remove l from Lρ if |tparent(l)|−|tl| ≤Mdiff,otherwise remove tl from tρ. Mdiff is a user-defined constant that intuitively results in merging cladesthat are too similar. At the end of the tree traversal, the set Lρ contains new candidate roots foreach initial clade. For each l ∈ Lρ define the summary CN profile as a vector whose i-th elementis the median of the copy number states of the i-th bin for all cells in tl. Compute the distancebetween two subclades as the mean absolute difference of their median genotypes. Merge subclonesinduced by Lρ if their summary CN profiles are too similar. We can do this by computing a t-testover the pairwise distances to exclude outlier subclades and merge the rest.Once clones are identified, we set the abundance of each clone at a specific timepoint as thefraction of cells in that clone from that timepoint. We note that for the data from WGS bulksequencing [43] we used the following procedure to estimate clonal fractions: (i) let ν denote themutational cellular prevalence (rows) estimated over multiple timepoints (columns) using the multi-sample PyClone [22] model, (ii) define β as the genotype matrix (which mutation-cluster (rows) ispresent in which clones (columns)), (iii) then we set βγ = ν where γ = β−1ν are the clonal fractionsover time, and (iv) we solve for γ using QR-decomposition.2.7 ConclusionsIn this chapter we introduced the sitka analysis pipeline, and described a number of visuali-sation methods that supplement the sitka phylogenetic reconstruction model. We briefly re-42viewed parts of the sitka model’s mathematical formulation that was related to the visualisationtoolkit. We also compared sitka over three real-world datasets to some benchmark methods usinga posterior-predictive test that we derived. We then introduced a tree-cutting procedure that wecall Lumberjack that is used to identify major clades in the trees inferred by sitka.In this work we use data in which the genome of the single cells CNA profiles are partitionedinto fixed-500Kb bins. Each bin is assigned a constant integer CN state. The size of these binsindicates that we cannot preclude that the same bin harbours multiple CNA events. Biologicalprocesses that result in complex DNA rearrangements could further increase the probability of thispossibly [69, 70]. These cases will violate the perfect phylogeny assumptions.It may be possible to take a data-driven approach to determining bin sizes. One approach wouldbe to use smaller bin sizes for regions of the genome where more sequenced reads are available.In such regions, the copy number states can be more accurately determined for the smaller bins.While in this work we have used the CNA states of pre-defined fixed size bins as input, segmentingthe genome, and calling the copy number jointly with inferring the phylogeny may be possible,although at a higher computational cost.Characterising existing subpopulations in a population of cells is a major goal of this work. Thepreprocessing step removes multiple cells (up to 90% of the sequenced cells, see Table A.5). Wefilter out a fraction of cells to remove contaminated cells, either doublets or mouse cells, cells withtoo many erroneous sequencing artefacts, and cycling cells. Removing a portion of the sequencedcells will decrease the statistical power to determine the subclonal structure of the population. Ourcurrent cell filtering step is focused on removing cycling cells (Section A.4). This may bias oursampling against clones that have a higher division rate. It would be beneficial to retain moresequenced cells without negatively impacting the quality of the downstream analysis. One possibleway is to retain false positive doublet cells by reviewing the images taken by the DLP+ platform’smicroscope prior to sequencing. As each library may contain up to thousands of cells, this approachis labour intensive and will benefit from automatic image processing frameworks. Another area toexplore is to decrease the error rate in the cell quality classifier [20]. This may be done by taking adata-driven approach to choosing the features used for scoring cells in the classifier and to retrainthe classifier on a larger dataset. A related idea consists in quantifying the error rate in differentregions of the genome; if in a cell with a low quality score, the genomic bins can be classified43into high-quality and low-quality regions, it may be possible to use the high-quality regions asphylogenetic markers and assign the cells to an already inferred phylogenetic tree.Structural variations such as Chromothripsis that affect multiple segments of the genome atthe same time make it difficult to determine the rate of CNA events and suggest that CNA eventsmay not be suitable molecular clocks to estimate branch lengths. One possible remedy is to firstinfer the tree topology via markers based on CNA events and then conditioned on this topology,add SNVs to the tree. The number of SNVs on each edge of the tree may give an indication of themutation rate in each clade.We note that the trees inferred in Figure 2.12-a-c are unbalanced. Unbalanced tree topologieshave been observed and are expected in adapting populations [71]. However this may be partiallydue to the fact that the copy number states of each cell is determined individually and independentlyfrom the tree inference step and that our jitter fix heuristic method cannot completely correct theCN calls. It may be improved by the joint inference of the CN states, the genomic segments, andthe tree topology.The visualisation toolkit developed here is static. Adding new annotations requires re-renderingthe entire plot. It would be useful to create an interactive toolkit where individual nodes on the treecould be selected to view more detailed annotations. This may include relevant genomic informationsuch as genes affected by copy number changes and if available, data on point mutations.Evaluating the performance of a phylogenetic reconstruction method over real-world datasetsis difficult, mainly due to a lack of ground through. One promising area of research is the useof CRISPR-Cas9 based lineage tracing [72]. In absence of ground truth data, we developed apredictive test (Section 2.5.1) that to our knowledge enables a first of a kind benchmarking ofphylogenetic inference methods over real-world SCG CNA datasets.Phylogenetic tree reconstruction is a principled way to identify subpopulations in a heteroge-neous single-cell population. In Chapter 4, we use sitka to infer phylogenetic trees in multiplemodel systems from which scDNAseq is available through the DLP+ platform. We then use atree-cutting algorithm to cut the sitka trees and assign single-cells to clones. This enables usto track the abundance of subpopulations over multiple timepoints. In the next chapter we in-troduce fitClone, a statistical framework that uses the timeseries clonal abundance data to inferevolutionary parameters of clones in model systems.44Chapter 3Modelling Fitness in LongitudinalData Via an Approximation to theWright-Fisher Diffusion Model3.1 IntroductionAn important but open question in the treatment of patients with cancer is why and how resistanceto therapy is developed. In many patients, despite initial positive response to therapy, the cancereventually relapses. Advances in next generation SCS as well as longitudinal xenoengraftmenthas enabled more accurate, quantitative measurements of tumours. By serially transplanting atumour sample into highly immunodeficient mice, a system called patient derived xenograft, onecan continuously monitor how tumour composition may evolve in a patient over time. This processgenerates timestamped samples comprising thousands of cells that through SCS can be measuredat thousands of genomic features including CNAs. This is in contrast to the standard of care whereoften only up to two timepoints are acquired and measured at an aggregate course level. Thenew higher resolution datatype introduces major challenges in extracting meaningful knowledgefrom the data, including (i) how to identify biologically meaningful groups of cells (i.e., clones)across multiple timepoints, and (ii) how to quantitatively reason about the underlying evolutionaryforces acting on the clones via their observed dynamics. In Chapter 2, we developed sitka andLumberjack to address the first question. Here we set out to tackle the second question.Tracking the relative abundances of the clones over time produces a timeseries reminiscent ofallele frequencies. Over time evolutionary forces act on a population, changing the allele frequenciesthrough adaptation, mutation and neutral evolution [19, 73]. In particular, clonal abundances are45similar to the fraction of the individuals in a population that harbour a specific allele. We adoptedideas from the population genetics literature to model the evolution of tumour growth via a statespace model, the dynamics of which is dictated by the diffusion approximation to the Wright-Fisher process. The algorithm infers evolutionary fitness coefficients for each clone by simulatingtrajectories consistent with their observed clonal abundances. To the best of our knowledge, theonly exact algorithm for the multivariate case (where we track more than one clone at a time)has infinite expected run time [74]. Therefore we decided to develop an approximation and itsrespective inference engine. These models will quantify whether observed dynamics are significantlydetermined by selection versus purely driven by random stochastic processes (i.e., genetic drift) -a key determinant on the path to predictive models. Those clones under selection then providehypotheses over mechanisms of disease progression through analysis of mutational or molecularfeatures.We describe in this chapter, two approaches to fitness modelling. The first is a Bayesian state-space model (fitClone) based on the Wright-Fisher diffusion with selection (Section 3.1.1). Thesecond is a deterministic logistic growth model (Section 3.9). We compare the two methods inSection 3.10. The comparison favours the Bayesian model, hence this is the model used in ourresults unless we specify otherwise.3.1.1 fitClone: a Bayesian fitness model for timeseries dataWe developed a Bayesian model and associated inference algorithm based on a diffusion approxi-mation to K-allele Wright-Fisher model with selection. We start with timeseries clonal abundancemeasurements over a fixed number of clones and estimate two key unknown parameters of interest:fitness coefficients si for clone i which represents a quantitative measure of the growth potential ofa given clone; and distributions over continuous-time trajectories, a latent (unobserved) populationstructure trajectory in ‘generational’ time.After briefly reviewing and setting notation for Wright-Fisher diffusions with selection (Sec-tion 3.2), we introduce the Bayesian model we used to infer quantitative fitness of clones fromtimeseries data (Section 3.3). We then describe a novel algorithmic extension that scales simul-taneous inference to dozens of clones (Section 3.5), ancillary methods for effective population sizeestimation (Section 3.6), and reference clone selection (Section 3.8).46Statistical methods for inference in the Wright-Fisher model is reviewed in [75]. In the bi-alleliccase, one method is to estimate Ne first and then use approximate Bayesian computation [76] totarget the posterior of S, using trajectories sampled from a Binomial distribution [77]. It is possibleto use the continuous transformation of the discrete Wright-Fisher model in which both time andstates (number of alleles/individuals) are scaled by a multiple of Ne, and then approximate it usingdiscretisation in both time and allele-frequency states [78–80]. An interesting direction is the exactsimulation of diffusions [81, 82]. For the one-dimensional Wright-Fisher diffusion with a generaldrift function, an exact simulation algorithm exists, although each sampled trajectory can take upto minutes to generate [83].3.2 Wright-Fisher diffusions with selectionWe take the convention of using uppercase symbols to denote random variables and lowercasesymbols for realisations of random variables. Let K denote the number of clones and denote byZt =(Z1t , . . . , ZKt)the relative abundance of each of the K clones at time t in the population.The process Zt satisfies, for all t, the constraints∑Ki=1 zit = 1 and zit ≥ 0 for i ∈ {1, . . . ,K}. Inother words it is a (K − 1) simplex. We model the process Zt using a Wright-Fisher diffusion withselection.A Wright-Fisher diffusion can be written in stochastic calculus notation asdZt = µs,Ne(Zt)dt+ σ(Zt)dWt (3.1)where {Wt} is a K-dimensional Brownian motion, and the functions µ and σ, defined below,respectively control the deterministic and stochastic aspects of the dynamics where we assume thismaps to the selection and genetic-drift components of the observed data. For z =(z1, z2, . . . , zK),the vector-valued function µs,Ne : RK → RK is defined asµs,Ne(z) =(µs,Ne1 (z), . . . , µs,NeK (z))µs,Nei (z) = Nezi(si − 〈s, z〉),where 〈x, y〉 is the inner product of vectors x and y; Ne, the effective population size, discussed in47more details in Section 3.6; and the parameters s = (s1, s2, . . . , sK) are called fitness coefficients.The interpretation of the fitness parameters is that if si > sj , then subpopulation i has highergrowth potential compared to subpopulation j. In practice, small values of si − sj can result in alarge difference in the observed clonal abundances (e.g., see clones A and B from the HER2+ PDXtimeseries analysed in Chapter 4). The matrix-valued function σ : RK → RK×K is defined asσ2(z) = [σ2i,j(z)]i,j∈{1,...,K}σ2i,j(z) = zi(δi,j − zj)where δi,j is the Kronecker delta. Given an initial value z, we denote the marginal distribution ofthe process at time t by Zt ∼WF(s,Ne, t, z). Equation (3.2) shows an example of a two dimensionalWF process with s = (s1, 0.0), Zt =(Z1t , Z2t), and Wt =(W 1t ,W2t).dZ1t = NeZ1t−1(s1(1− Z1t−1))dt+√Z1t−1(1− Z1t−1)dW 1t (3.2)3.3 fitClone model specificationGiven timeseries data measuring the relative abundances of K subpopulations at a finite number oftimepoints, the output of the fitClone model is a posterior distribution over the unknown parame-ters of interest: the fitness parameters s described in the previous section, and the continuous-timetrajectories interpolating and extrapolating the discrete set of observations.To do this, fitClone places a prior on the fitness parameters s, uses a state space model (Fig-ure 3.1) in which the latent Markov chain is distributed according to a Wright-Fisher diffusion, andencodes the noisy sampling from the population at a discrete set of timepoints via the observationmodel.48Z0 Z1 Z2 ZTY1 Y2 YTθFigure 3.1: A state space graphical model. The observed quantity Yt denotes the noisy clonalabundances while Zt denotes the latent clonal fractions. θ = (Ne, S) comprises the parameters ofthe model. The Wright-Fisher process controls the transition between states.We model the fitness parameter with a uniform prior over range I:Sk ∼ Uniform(I), k > 1,where we set S1 = 0 to make the model identifiable (see Section 3.8 for details). We used I =(−10, 10) in our experiments. Note that the posterior is contained far from the boundaries of thisprior range in all experiments.The initial distribution, i.e., the distribution of the value of the process at time zero, is aDirichlet distribution with hyper-parameter (1, 1, . . . , 1),Z0 ∼ Dirichlet(1, 1, . . . , 1).This can equivalently be seen as a uniform distribution over the K-simplex.Let t1 < t2 < · · · < tT−1 < tT denote a set of process times at which measurements of relativeabundance of clones are available. Ideally, we would like the latent transition kernels to be given49by the marginal transitions of the Wright-Fisher diffusion from the last section,Ztm |Ztm−1 , S ∼WF(S,Ne, tm − tm−1, Ztm−1), (3.3)where Ne is estimated as a preprocessing step (Section 3.6). In practice we resort to approximatingthe distribution in Equation (3.3) via an Euler-Maruyama scheme (Section 3.4.1).Finally, for each t ∈ {t1, t2, . . . , tT }, let Yt =(Y 1t , . . . YKt)denote a noisy observation of theclonal prevalences at process time t. In the single-cell context, this is obtained by counting, for eachclone defined in Section 2.6, the number of cells coming from each passage, and normalizing by thenumber of cells sequenced in that passage. In the bulk sequencing context, see Section 2.6. Whenobservations are counts, a Multinomial observation model or the Normal distribution approximationto it can be used. For simplicity, in both cases we use a Normal observation model, i.e., yit | zit ∼N (zit, σ2obs), where σ2obs = npi (1− pi) and n = ∑j yjt and pi = yit/n.3.4 Posterior inference under the fitClone modelSince the marginal distributions of the Wright-Fisher diffusion do not admit closed form expressions,and previous work on exact simulation does not scale to high values of K, we resort to discretisationusing an Euler-Maruyama scheme [84].3.4.1 Euler-Maruyama discretisationIn datasets available to us, the latent state of the Wright-Fisher process, i.e., the relative clonalfractions are observed at a few timepoints. In other words the system is sparsely observed. Notethat the Wright-Fisher diffusion is a continuous time process where time is measured in units ofNe. Time represents the number of generations scaled by 1/Ne and each generation comprises acell division event. However, in order to perform inference using this model we need to discretisetime. Intuitively, this discretisation in time can be thought of as an intermediate approximationbetween sparse observations.We discretise the interval (tm, tm+1) usingN equidistant intermediate timepoints(tm,1, tm,2, . . . , tm,N)50such that h = tm,i − tm,i−1 = (tm+1 − tm) /N and h is the discretisation step size. For simplic-ity, we used the same h between all observed time steps and moved the observed timepoints totheir closest grid point. We augment the state space of the process Zt and add between Ztm andZtm+1 , N imputed intermediate states(Ztm,1 , Ztm,2 . . . , Ztm,N). Using these intermediate states,the marginal transition kernels in Equation (3.3) can be written as:p(Ztm+1 |Ztm , S)=∫. . .∫p(Ztm+1 |Ztm,n , S) { N∏i=2p(Ztm,i|Ztm,i−1, S)}p(Ztm,1|Ztm , S)dZtm,1 . . . dZtm,Nwhere the intermediate transition kernels, i.e., the Euler-Maruyama steps are approximated viaa Normal distribution as follows:Ztm,i |Ztm,i−1 , S ≈ N(Ztm,i−1 + µS,Ne(Ztm,i−1)h, σ2(Ztm,i−1)h)(3.4)In Chapter 4 we apply this discretisation to timeseries data from in vitro and in vivo systems.For instance, in Section 4.3, a TNBC PDX timeseries dataset passaged over 1,002 days is analysed.We estimate the diffusion time t10 − t1 = 0.90 and use a discretisation step of h = 0.05 to yield 8intermediate steps.3.4.2 InferenceWe used a particle Markov chain Monte Carlo (pMCMC) method called Particle Gibbs with An-cestor Sampling and particle rejuvenation [85, 86] to sample from Ztm and the intermediate Euler-Maruyama steps, and a Metropolis within Gibbs sampler to sample the selection parameters S.In sampling S, and for a fixed trajectory z = (zt1 , zt2 , . . . , zT−1, zT ), we use a truncated Normalrandom walk proposal with variance σp and boundaries (Il, Iu) that are set to the boundaries ofthe prior on S. This proposal is non-symmetric and yields the following MH acceptance ratio:a = min[1.0,p(s′|z)ΦT N (s|s′, σp, Il, Iu)p(s|z)ΦT N (s′|s, σp, Il, Iu)]where s′ denotes the proposed value and ΦT N(. | s, σp, Il, Iu)is the probability density function51of a truncated Normal distribution parameterised by its mean, variance, and lower and upperbounds.We note that from p(s|z) = p (z|s) p(s)/p(z) only the likelihood term p (z|s) remains aftercancelling identical factors in the numerator and denominator of the acceptance ratio above. Thislikelihood term factors into marginal distributions of the Wright-Fisher diffusion:p(z|s) = T−1∏t=1p(zt+1|zt, s)In practice we keep the Euler-Maruyama intermediate steps and Equation (3.4) can be used toevaluate p(z; s) for a given s.3.5 Conditional samplerIn sampling a trajectory between a pair of consecutive latent states ztm−1 and ztm , the samplerinitialises a set of particles and evolves them according to the distribution of the Euler-Maruyamaintermediate steps to form potential trajectories. Each particle represents a potential latent state.A trajectory is valid only if it ends inside a K-dimensional sphere centred at the observation attime tm. We call this sphere an balltm . This is essentially a soft-bridging problem and becomes moredifficult with increasing K as more components of a potential trajectory have to pass through theballtm . Given a valid trajectory between ztm−1 and ztm , it is possible to evolve only a few componentsof the potential trajectory conditioned on the other components remaining constant. In other words,instead of sampling the clonal fraction of all clones at the same time, we keep the abundance ofa subset of clones fixed and conditioning on those, update the clonal fraction of the rest. Thediscretised Wright-Fisher diffusion transition kernel at time t is normally distributed so we can usethe conditional distribution form of a multivariate Normal distribution.Without loss of generality we focus on a single intermediate Euler-Maruyama transition step,namely Ztm,i |Ztm,i−1 , S. In the following it is useful to distinguish the subset of parameters that aresampled from those that are fixed. Let B = (B1, B2) partition the index set 1:K = {1, 2, . . . ,K}into two subsets, that is B1 ⊂ 1:K, B1 ∪ B2 = 1:K and |B1| 6= 0. We denote by B1 the indexes ofthe subset of parameters that are being sampled.52Define ZBjtm:={Zitm | i ∈ Bj}, then we can write Ztm =ZB1tmZB2tm. Now the problem can bestated as sampling from the following conditional distribution:ZB1tm+1 |ZB1tm , ZB2tm+1 , ZB2tm , SWe use a two step procedure to achieve this: First, we sample Z¯B1tm+1 fromZ¯B1tm+1 |ZB1tm , ZB2tm+1 , ZB2tm , SSecond we ensure ZB1tm+1 respects the simplex criterion and set a realisation of ZB1tm+1tozB1tm+1 = fzB2tm+1(Z¯B1tm+1)where the vector-valued function fy : [0, 1]J → [0, 1]J for the vectors x and y of sizes J and K − Jis defined as fy(x) =(f1(x; y), f2(x; y), . . . , fJ(x; y))and each component is equal tofj(x; y) =∣∣xj∣∣ , j < i,max[0, 1−∑K−Jr=1 yr −∑i−1r=1|xr|] , j = i,0.0 otherwise.wherei = argmini∈{1,...,J} i∑r=1|xr| ≥ 1.0−K−J∑r=1yr .The effect of the mapping f(.) becomes negligible as the discretisation step size goes to zero,but it is needed in a finite discretisation to ensure that the simplex constraint is satisfied. Nowwe describe sampling from Z¯B1tm+1 . To reduce clutter we drop the dependence on parameters s andNe and assume they are fixed to some given value. Then, from the conditional distribution of amultivariate Normal we get:Z¯B1tm+1 | ZB1tm , ZB2tm+1 , ZB2tm , S ∼ N(µ¯, Σ¯)53whereµ¯ = µB1 + ΣB1,2(ΣB2,2)−1 (zB2tm+1 − µB2)andΣ¯ = ΣB1,1 − ΣB1,2(ΣB2,2)−1ΣB2,1Σ¯ is the Schur’s complement of ΣB2,2 in Σ.µ (ztm) = µ =µB1µB2 (3.5)Σ (ztm) = Σ =ΣB1,1 ΣB1,2ΣB2,1 ΣB2,2 (3.6)where in a slight abuse of notation, ΣBi,j :=[Σp,q]p∈Bi,q∈Bj . For a more efficient computation,we precompute and cache the value of(ΣB2,2)−1the inverse of ΣB2,2 . For data analysed in thiswork, we only update one component at a time, that is we set |B1| = 1.3.6 Estimating the effective population sizeFollowing [87] we use F′s an unbiased moment-based estimator of the Ne where Ne =1F′Sand t isthe number of generations between each passage.F′s = (1/t)Fs(1− 1/(2n˜))− 1/n˜(1 + FS/4)(1− 1/ny) (3.7)where Fs =(x−y)2z(1−z) and z = (x + y)/2 and n˜ =2nynxny+nx, the harmonic mean of the sample size(initial population size at the passage) nx and ny at the two timepoints. x and y are the minorallele frequencies at the two timepoints.In the multi-allelic case, we have:Fs =1KK∑i=1(xi − yi)2zi (1− zi)54This is equivalent to plan 2 in [87], sampling before reproduction and without replacement. Inour formulation, the multi-allelic maps to multi-clones and this is what we have used to estimateNe in this work.We used the sum of clone sizes as the approximate initial population size at each timepoint/passage.Table A.3 lists the resulting Ne estimates. Since fitClone is robust to the choice of Ne in this range(Figure 3.2-c), we set Ne = 500 for all datasets analysed in this work. A smaller Ne results in morestochasticity in the estimation of selection coefficients. We have chosen this value for Ne in partto account for phenomena that we do not explicitly model in our formulation of the Wright-Fisherdiffusion, including clonal interference. We do not claim that Ne = 500 represents the populationstructure.We note that in our model we assume that the effective population size remains constant overall timepoints. This does not take into account the potential changing population growth rate orthe bottleneck effect due to passaging. These phenomena may scale the diffusion time and biasour estimates of evolutionary events, including fixation or extinction times (see also Section 4.3.2).This stretching and compressing of time could be accounted for by adding random effects to thenumber of generations in the model, for example by taking Ne,t as a piece-wise constant randomvariable that can vary between passages. Informally, the estimator introduced in this section canbe thought of as a corrected harmonic mean estimator of the census population sizes over timethat can to some degree account for the chaning population size [88]. Table A.4 shows the otherparameters used in the inference over the real-world datasets.3.7 Summarising the posterior distributionIn all real-world data that we analysed in Chapter 4, 10,000 particles and a burn-in equal to 10%of the MCMC samples were used. For the trajectories, we reported zˆ1:T where zˆt =(zˆ1t , . . . , zˆKt)and zˆkt encodes the posterior mean of the clonal fraction of clone k at time t.The posterior of the selection coefficient vector was summarised by sˆ = (sˆ2, . . . , sˆK) where sˆkdenotes the posterior mean of the selection coefficient of clone k (see Figure 4.9). To comparethe selection coefficients of two clones we used a (K − 1) × (K − 1) posterior ordering matrixP (Figure 3.2-b,c). Pi,j = P (Si ≤ Sj | y) shows the posterior probability that clone i has a55higher selection coefficient than clone j. We estimate this quantity by Pi,j =∑Mm=1 1(sm,i > sm,j)where sm,k is the sampled selection coefficient of clone k at MCMC iteration m. Figure A.6 showsthe visualisations of the posterior ordering matrices (P ) for the real-world datasets analysed inChapter 4. In these visualisations we colour coded the value of Pi,j , with the stronger purplehues (close to 1.0) representing a higher confidence that clone i dominates clone j, and converselythe stronger grey hues (close to 0.0) denote that clone j dominates clone i. Colours closer towhite (0.5) represent no dominance. Note that for the lower diagonal elements P(Sj ≤ Si | y)=1− P (Si ≤ Sj | y) and are omitted for clarity. The diagonal entries are to guide the eyes only.3.8 Selecting the reference cloneIn our formulation of the Wright-Fisher diffusion, one reference clone with a selection coefficientof zero has to be chosen. The selection coefficient of the other clones are reported relative to thisvalue. For instance, if the fittest clone is chosen as the reference, the other clones will have negativeselection coefficients. We chose to set the reference to a clone with an approximately monotonicallydecreasing trajectory (clonal abundance over time). This choice was motivated by a desire to infera non-negative value for the fittest clones. Figure 3.2-b shows that the model is robust to thechoice of the reference clone. We ran the inference procedure over the same dataset multiple times,each time changing the reference. The posterior ordering of clones over different choices of clonesremained mostly identical.3.9 The deterministic logistic growth modelWe developed a closely related population genetics model which incorporates selection via deter-ministic differential equations, but has closed form solutions. In this model [36], the solution of adeterministic DE, the frequency of each population c out of possible K populations, at time t, withselection coefficient sc is proportional to its starting prevalence p0(c) multiplied by a power of therelative fitness coefficient wc = (sc + 1),f(c, t, s) =(sc + 1)tp0(c)∑Kk=1(sk + 1)tp0(k)(3.8)56To estimate fitness coefficients w1:T from observed clonal fractions y1:T , we solve the optimizationproblem in Equation (3.9) using a limited memory Broyden-Fletcher-Goldfarb-Shanno (BFGS)optimization procedure with box constraints [89]. Note that ft,c(w) is the clonal fraction estimatefrom Equation (3.8) and yt,c is the observed clonal fraction for clone c at time t.minw∈W T∑t=1dt(w) (3.9)where d(w) =(d1(w), d2(w), . . . , dT (w))is a vector valued function whose elements are dt(w) =(∑Kc=1(ft,c(w)− yt,c)) 12and W = (R∪ {0})K and w1 = 1.0.3.10 Simulation benchmarkingWe forward simulated L = 40 datasets from the joint distribution of the Wright-Fisher modelwith K = 5, and L = 40 datasets with K = 11. For each simulated dataset l, we sampled theinitial clonal abundance vector Z1,l ∼ Dirichlet (α1:K) where αi = 1 and selection coefficients froma Normal distribution with si,l ∼ N (0.0, 0.3) truncated at (−0.5, 1) for i ∈ {2, . . . ,K} assumingthe index of the reference clone is i = 1 and s1 = 0. Discretisation constant (step size) ∆τ =0.001 and the standard deviation of the emission model was set to σobs,simul = 0.001. At simulationand inference, we set Ne = 500. The simulation was continued to a diffusion time of 0.1 afterwhich 10 equidistant samples were recorded as observed values for the process. In all modelsexcept the Logistic growth, we put a uniform prior on each component of the s vector, that is,si ∼ Uniform(−5.0, 5.0), and a Dirichlet prior on the initial clonal distributions. We set step size∆τ = 0.001, σobs,infer = 0.01 and used 10,000 particles for 10,000 MCMC iterations.We ran 5 different models on the simulated dataset as follows (Figure 3.2): (i) WFda is theWright-Fisher model with diffusion approximation. (ii) Logistic growth is the deterministic differ-ential equation Wright-Fisher model. (iii) One-step is identical to the WFda model but applies nodiscretisation for the time between observations. (iv) Single Clone is the WFda model with K = 2.See Section 3.9 for how the logistic model was fit. For the other models, we used the mean absoluteerror (MAE) of the post burn-in mean posterior of marginal selection coefficients for each methodand the selection coefficient used to generate the simulated data. For each dataset l, the Single57Clone model is run K times, once for each clone k, where its input consists in the observed clonalfraction of only one clone. For this model we reported the averaged MAE across clones per dataset.WFdaWFda5 11acbFigure 3.2: Simulation studies for the fitClone model. (a) Comparison to baseline methods forK = 5 clones (left) and K = 11 clones (right). (b) Posterior ordering of clones based on theirinferred posterior selection coefficients across three values of effective population size (columns)and different choice of the reference clone (rows). (c) Posterior ordering of clones across differenthyper parameters in the fitClone model. Effective population size in the range of 500 to 5,000(top column), and observation error (bottom column). Minimum number of interpolations betweentwo observations (rows).The WFda and the logistic growth methods result in comparable error profiles, with a slightadvantage for the WFda method in the K = 11 regime. This may be due to accounting forstochasticity in the dynamics via the diffusion process. Since the logistic growth is our baseline, aperformance in line with this method suggests that our implementation of WFda is correct.Recall that in the Euler-Maruyama scheme, N denotes the number of discretisation steps be-58tween two consecutive timepoints (tm, tm+1) where clonal abundance measurements are available.The One-step method uses N = 0. When K is small, this setting is in line with that of WFda andlogistic growth. This suggests that with fewer clones, a coarser discretisation (smaller N) can beused. However, when K is large, the One-step setting results in an error distribution with a largevariance where the inferred selection coefficients in a large subset of datasets have higher meanabsolute errors. This suggests that even though the computational cost of inference rises in generalwith increasing the number of clones, decreasing the discretisation budget may result in inaccurates estimates and should be avoided.The Single Clone model has the highest estimation error. This suggests that the K = 2 andK > 2 regimes can generate similar clonal abundance trajectories with different underlying s values.Therefore treating each clone independently when estimating s is sub-optimal.Comparing K = 5 to K = 11 across the methods compared here, we note that increasing thenumber of clones shifts the error distribution for s upwards. This may be due to unidentifiabilitywhere as K increases, more configurations of s are compatible with the observed clonal abundancetrajectories. Observing the same timeseries for longer, or having more parallel repeats of thetimeseries, for instance, establishing a parallel branch, may help improve the estimation accuracy(see Chapter 4).Figure 3.2-c shows the model fits over a range of parameters for a dataset with K = 5. Werepeated the inference procedure with Ne ∈ {500, 1000, 1500, 2000, 2500, 3000, 5000}, and the obser-vation error σ2obs ∈ {0.025, 0.05}, andN∗ ∈ {1, 4} whereN∗ is the minimum number of discretisationsteps between two consecutive timepoints. The posterior ordering matrices show that in generalthe model is more confident in the difference between the s coefficients as the Ne rises. Overallthe model fits were robust to the choice of the effective population size and the discretisation stepsize. Together these simulations established a rationale for systematic modelling of all clones in aunified approach with a generative process.3.11 Proof of principle application on TNBC-SA501 PDXWe applied fitClone to previously published data, wherein experimentally derived reproducibleclonal dynamics had been reported in breast cancer PDXs [43]. From one of these lines, the59abundance of five major clones (A-E), as determined by single cell genotyping was measured oversix serial passages. Note that in this dataset, unlike those analysed in Chapter 4, clones wereidentified by their point mutation profiles. One clone (clone E) was found in the original studyto undergo a selective sweep in repeat and independent passaging [43], suggesting selection due tohigher fitness.The input to fitClone are the relative abundances of the five clones at every timepoint, alongwith an estimate of the number of generations up to each timepoint (Figure 3.3-a). From this,we infer Ne as described in Section 3.6. For each timepoint, we set the number of generations tothe number of days from tumour transplantation to tumour collection for in vivo systems, and thenumber days in culture for in vitro systems. In Section 2.6 we explain how the clonal fractions werederived from the bulk WGS data in more detail.fitClone projects the observations from the generation time into the diffusion time which is inunits of Ne (Section 3.4.1). It then applies the Euler-Maruyama scheme to discretise time and im-putes the clonal fractions at unobserved timepoints between the observed timepoints. The outputsof the model include the clonal fractions at imputed timepoints consistent with the observations,and their estimated selection coefficients. Figure 3.3-b shows the observed and imputed clonalfractions in diffusion time. The dashed white vertical lines correspond to the discretised timesat which clonal fractions were observed. Figure 3.3-c,d,e shows the distribution of the estimatedselection coefficient of each clone, the credible intervals, and the posterior ordering matrix.fitClone estimates converged with Clone E bearing the highest fitness (1 + s=1.03 ± 0.01),consistent with positive selection over the timeseries (Figure 3.3), and thus representing a proof ofprinciple application of fitClone on real-world data.60Figure 3.3: The fitness analysis for clones derived from bulk WG. The dataset is SA501. Applicationof fitClone to previously published clonal dynamics in PDX showing observed data (a), inferredtrajectories (b), posterior distributions of fitness coefficients (c), the credible intervals for fitnesscoefficients (d) and pairwise comparison of the selection coefficients (e).3.11.1 Software and implementationThe software implementation of fitClone is available at: [https://github.com/UBC-Stat-ML/fitclone]3.12 ConclusionsIn this chapter, we introduced a statistical framework, fitClone, that models timeseries clonalabundance observations using an implementation of the Wright-Fisher diffusion process. It si-multaneously estimates growth trajectories, Zi and fitness coefficients, Si for each clone i in thepopulation (Figure 3.3). We note that increasing 1 + Si indicates positive selection and highergrowth potential. The model accounts for drift as well as selection, with fitness estimated relativeto a reference population, where S = 0 by construction. As a generative process, the model can beused for forecasting evolutionary trajectories of specific clones. In the next chapter, we will applyour framework to multiple model systems of cancer, measured at the single cell level. We will use61sitka to infer the clonal structure of heterogeneous single cell populations, and then will applyfitClone to model the clonal dynamics.f62Chapter 4Application of the Wright-FisherDiffusion Approximation to CancerModel Systems Interrogated at theSingle Cell Level4.1 IntroductionTumour fitness landscapes underpin selection in cancer evolution and response to treatment. How-ever, quantifying fitness in heterogeneous cell populations remains an open problem that hindersprogress in developing effective therapeutic strategies. In Chapter 1, we posited that two key tech-nological advances, namely the emerging single cell sequencing platforms and model systems, inparticular longitudinal patient derived xenoengraftment, may in principle facilitate the monitoringof change in tumour cellular composition over time and in response to therapeutic interventions.Measurements taken from these systems constitute novel datatypes that introduce major analyticalchallenges including (i) how to identify biologically meaningful groups of cells (i.e., clones) acrossmultiple timepoints, and (ii) how to quantitatively reason about the underlying evolutionary forcesacting on the clones via their observed dynamics. In Chapters 2 and 3, we described two main com-putational ingredients to address challenges above, namely sitka, a phylogenetic inference methodand accompanying toolkits to assign cells to clones, and fitClone, a probabilistic framework thatascribes quantitative selection coefficients to individual cancer clones and forecasts competitiveclonal dynamics over timeIn this chapter we apply the methods developed in the previous chapters to real-world datasets63to make inferences about how human cancer cells evolve at the copy number level, and that howestablishing baseline fitness measures helps to interpret selection under drug administration. Wehave developed an experimental and computational platform consisting of three major components:scaleable phylogenetics for single cell genomes to identify clones (sitka and Lumberjack), timeseriessampling of immortal cell lines and patient derived xenografts to observe clonal dynamics, and amathematical model for inferring clone-specific fitness measures (fitClone). Figure 4.1 shows theexperimental setup and the data analysis, from the experimental setup in the cell lines and the PDXtimeseries up to the clonal fraction measurements over time, the input to the fitClone model.For normal human breast epithelial cells [90] in vitro and in breast cancer PDX [43, 91] (Fig-ure 4.1-a), we sequenced >60,000 cells over interval passaging (Figure 4.1-b) with single cell wholegenome sequencing [20] measuring single cell copy number profiles, and computing phylogenetictrees (via sitka) to identify genotypic clones (using Lumberjack) and their relative abundances asa function of time. Genetic (p53 biallelic inactivation for cell lines) and pharmacologic (cisplatindosing in PDX models) perturbations were applied to determine their impact on fitness landscapes.fitClone was used to measure the selection coefficient for each clone (s) which we hypothesise toindicate growth potential. The larger the value of s, the more fit the clone is relative to the chosenreference clone.The chapter is organised as follows. We first investigate the evolution of clonal trajectories invitro in three related cell lines (Section 4.2) and in vivo with two independent breast cancer PDXtimeseries (Section 4.3.1). Next we turn our attention to in vivo experiments where three differentindependent PDX lines are analysed to (i) explore the predictive power of our model throughestablishing physical mixture lines (Section 4.3.2), and (ii) to determine the dynamics of selectionunder early response to cisplatin treatment (Section 4.4.2).64Figure 4.1: Schematic overview of experimental design for quantitatively modelling clone-specificfitness. (a) Timeseries sampling from in vitro or PDX systems. (b) Clonal dynamics of cellpopulations observed over time. (c) Whole genome single cell sequencing of timeseries samples.(d) Phylogenetic tree inference using sitka. (e) Clonal fractions over time constituting the inputto the fitClone framework, the mathematical modelling of fitness with diffusion approximationsto the K-type Wright-Fisher model.4.2 In vitro systems4.2.1 Serial passaging of immortalized 184hTert diploid breast epithelial celllinesWe applied our framework to immortalized 184hTert diploid breast epithelial cell lines [90] todetermine mechanisms by which TP53 mutation induces clonal expansions and fitness trajectories.Cell lines, relative to PDX systems, are easy to establish and perturb. In particular, they are wellsuited for high-throughput drug screening [92, 93]. We note that our analysis here is a proof ofprinciple and more replications are required to generalise our observations and conclusions. Weinvestigate three related timeseries (branches). Their relationship is shown in Figure 4.1-a.TP53 is the most abundantly mutated gene in all human cancers [94] and specifically in breast65cancers [10, 95]. Known to be permissive of genomic instability, TP53 loss is often acquired earlyin evolution and results in profound alteration of the copy number landscape [17, 96–98]. We askedwhether specific clonal expansions could be observed, and moreover if selective fitness advantagescould be quantified as a function of TP53 ablation, thereby modelling copy number driven etiologicprocesses in a controlled system of immortalized mammary epithelial cells. TP53 wildtype (p53WT ) timeseries sampling (60 passages over 300 days, 4 samples) was contrasted with two isogenicTP53 deficient NM_000546(TP53):c.[156delA];[156delA] [20] parallel branches (p53-/-a andp53-/-b), each passaged over 60 generations (285 and 220 days, respectively) and sampled 7 times.A median of 1,231 cells per passage were sequenced yielding a total of 6,620, 7,935, 9,615 single cellgenomes for each timeseries, respectively (Figure A.7-a,b).4.2.2 Phylogenetic and fitness analysisFor each of p53 WT, p53-/-a, p53-/-b we inferred single cell copy number profiles, constructed aphylogenetic tree to establish clonal lineages (Section 2.6, A.8) and measured clonal abundancesas a function of time (Figure 4.2-a,b,c). For each timeseries (branch), we pulled all single cellsand used sitka to infer a phylogenetic tree. The tree was cut using Lumberjack to yield clonesand their abundances over time. In the following, we summarise the fitClone inferred selectioncoefficient of a clone by reporting the posterior mean 1 + s followed by its standard deviation. SeeFigure A.5 for highest posterior density credible intervals for the selection coefficient of all datasetsthat we analyse in this chapter.Modelling the abundances with fitClone (Tables A.3, A.4) revealed p53 WT clonal trajectoriesconsistent with small differences over the posterior distributions of fitness coefficients amongst fourmajor clones (Figure 4.2-d). In contrast, p53 mutant branches each showed expansions of cloneswith aneuploid genotypes, where the founder diploid population was out-competed. Relative to p53WT, rates of expansion of p53 mutant, aneuploid clones were higher, leading to rapid depletion ofdiploid cells (Figure 4.2-d). Pairwise difference of selection coefficients si − sj (∆s) between clonesin the fitClone inference process was larger in the p53-/- lines relative to p53 WT (Figure 4.2-e).This suggests that p53 mutation permits the expansion of clones at higher rates, and that theseclones have measurably higher positive selection coefficients.The p53 mutant lines harboured 11 (size range 47 to 1,474 cells, median 204), and 10 (size66range 158 to 997 cells, median 404) distinct clones for p53-/-a and p53-/-b, respectively. In eachseries the diploid founder clones, devoid of detectable copy number alterations, were systematicallyout-competed by populations that had acquired at least one copy number alteration (Figure 4.2-panels f to m). Some similarities in copy number events were observed during the two replicatetimeseries. These included gains in chromosomes 13, 19p and 20, and losses on chromosomes 8p and19q (Figure A.8). However, by the end of the timeseries, the genotypes in the two lines had divergedconsiderably. Selection coefficients were highest in clones with localised amplifications of knownprototypic oncogenes in breast cancer [10, 95, 98, 99] including in MDM4, MYC (Figure 4.2-f) andTSHZ2 (Figure 4.2-j), in some cases on a whole genome doubled background. Clone A, the highestfitness clone in p53-/-a (57% of cells at last timepoint, 1+s = 1.05 ± 0.09) exhibited a whole genomedoubling event (18 chromosomes with four copies) and harboured a focal, high level amplificationat the MDM4 locus on Chr1q (Figure 4.2-f). Clone G (27% of cells at last timepoint, 1+s = 1.03 ±0.03), the next highest fitness clone in p53-/-a remained diploid, with the exception of a focal highlevel amplification precisely at the MYC locus on Chr8q (Figure 4.2-g). By contrast clone K, chosenhere as the reference clone for modelling (see Section 3.8), remained entirely diploid and exhibiteda monotonically decreasing trajectory (from 90% to 0% of cells over the timeseries, Figure 4.2-h,i).In p53-/-b, two clones exhibited non-neutral, positive selection coefficients (Figure 4.2-j). Clone D(52% of cells at last timepoint, 1 + s = 1.05 ± 0.02) harboured a Chr20q single copy gain with anadditional high level amplification at the TSHZ2 locus, while clone E (35% of cells at last timepoint,1 + s = 1.05 ± 0.04) harboured a Chr4 loss, Chr19p gain/19q loss and Chr20q single copy gain(Figure 4.2-f). As seen in p53-/-a, the ‘root’ clone I that remained diploid was systematicallyoutcompeted, diminishing from 68% to 0% abundance over the timeseries (Figure 4.2-k,l). Wenext tested whether the genomes of individual cells became progressively more aberrant over timeand whether this correlated with measurements of clonal fitness. We estimated both sample andclone specific mutation rates at the copy number breakpoint and point mutation level (as previouslydescribed [20]). Both p53-/-a and p53-/-b exhibited increased mutations and breakpoints over timerelative to p53 WT (Figure 4.3a-e). Cells accumulated 0.08 additional breakpoints (p<0.0001)and an average of 0.4 additional mutations (0.17 in p53-/-a and 0.67 in p53-/-b) per generation(Figure 4.3-k,l), while the p53 WT line accumulated 0.03 additional breakpoints per generation(Figure 4.3-k). This cell line was used as the reference for clone-specific point mutation detection67in p53-/- cell lines and as such we did not call or analyse any point mutations in p53 WT. Clonelevel distributions of breakpoints and mutations were positively correlated with inferred fitnesscoefficients in both p53-/- lines (Figure 4.3-g-j), but not in p53 WT (Figure 4.3-f).4.2.3 Segmental aneuploidies correlate with positive selection in diploid p53deficient cellsThese results indicate that the impact of genetic perturbation on selection can be measured andmodelled in clonal populations. In particular, p53 deletion, known to be an early event in theevolution of many cancers [100], yields clonal expansions driven by whole genome, chromosomaland segmental aneuploidies, conferring quantitative fitness advantages over cells that maintaindiploid genomes. Modes of positive selection involving high level amplification of proto-oncogenesoften seen in human breast cancer [101–104], and aneuploidies in general, suggest that in vitrogenetic manipulations can induce fitness-enhancing genomic copy number changes consistent withetiologic roles in cancer.68Figure 4.2: Impact of p53 mutation on modes of selection on 184hTert cells. (a-c) Clonal dynamics of p53 wildtype, and two independenttimeseries of p53 mutant 184hTert mammary epithelial cell lines. (d) Clonal fraction of diploid reference over time. (e) Distribution overmagnitude of difference between selection coefficients of pairs of clones. Each point is analogous to an element in the posterior orderingmatrix of Section 3.7. (f) Clonal genotypes of three representative clones for p53-/-a showing high level amplification of MDM4 in cloneA and MYC in clone G. Reference diploid clone K shown for comparison. (g) Phylogeny (simplified type II sitka tree) of cells over thetimeseries p53-/-a where nodes are groups of cells (scaled in size by number) with shared copy number genotype and edges representdistinct genomic copy number change points (sitka markers). (h) Inferred trajectories and (i) quantiles of the posterior distributionsover selection coefficients of fitClone model fits to p53-/-a. See Figure A.5 for the credible intervals and Figure A.6 for pairwisecomparison of the selection coefficients. (j) Clonal genotypes of three representative clones for p53-/-b showing high level amplificationof TSHZ2 in clone D, Chr4 loss in clone E. Reference diploid clone I shown for comparison. (k-m) Analogous to (f-i) but for p53-/-b.69Figure 4.3: Structural variant and mutation rates of 184hTert cells. Distribution over copy numberbreakpoints/cell as a function of generation for (a) p53 WT (b) p53-/-a (c) p53-/-b. Distributionover point mutations/cell as a function of generation for (d) p53-/-a and (e) p53-/-b. Clone specificdistributions over copy number breakpoints/cell, coloured by fitness coefficients for (f) p53 WT (g)p53-/-a and (h) p53-/-b. Clone specific distributions over point mutations/cell, coloured by fitnesscoefficients for (i) p53-/-a and (j) p53-/-b.4.3 In vivo system in the untreated regime4.3.1 Serial passaging of TP53 mutant PDX tumoursWe next modelled clonal expansions observed during serial passaging of TP53 mutant PDX tu-mours to determine the predictive capacity of fitness coefficients (Figure 4.4). PDX systems aremore difficult and time consuming to establish than cell lines, but may represent a more realisticmodel of tumour evolution as cells are located in their micro-environment. Here we exemplify ourframework on two PDX timeseries, namely HER2-positive-SA532 and TNBC-SA609, and contrasttheir evolutionary characteristics.We generated single cell genomes from 8 serial PDX transplants over 721 days from a HER2positive (HER2+) breast cancer with a TP53 p.A159P missense mutation (SA532), and contrastedthis with 10 serial samples over 1,002 days from a TNBC PDX (SA609) with a TP53 p.R213* non-sense mutation (Figure A.7-c,e). A median of 907 single cell genomes were sequenced per passage for70a total of 11,705 and 10,553 single cell genomes from the HER2+ and TNBC series, respectively.Both series exhibited progressively higher tumour growth rates over time (Figures 4.4-f,j, A.7-d). Data were analysed as per the in vitro lines described above and modelled with fitClone(Tables A.3, A.4). The HER2+ series exhibited 4 distinct clones ranging in size from 134 to 1,421cells (median 319, Figure 4.4-e), and the TNBC series exhibited 8 distinct clones with 18 to 680 cells(median 556, Figure 4.4-i). In the HER2+ model, clonal trajectories were consistent with selectioncoefficients with small relative differences in fitness (Figure 4.4-f,g). Clones B, C, and D showsimilar but distinct copy number profiles. The median genotypes of Clones B and C are differentacross 26 genomic bins, that of clones B and D are different in 24 genomic bins, and those of clonesC and D have differences along 21 genomic bins. By contrast, the TNBC model trajectories resultedin a high positive selection coefficient for a minority of clones (Figure 4.4-j,k). Consistent withincreased dynamics in the TNBC series, we found an initial increase of 0.1 breakpoints per cell pergeneration in the first 4 passages (Figure A.10-a). After this initial increase the average number ofbreakpoints per cell remained mostly constant. In the HER2+ line we observed a small decrease of0.04 copy number breakpoints per generation. We note that clone E in TNBC swept through thepopulation over the last 3 timepoints (Figure 4.4-j,k, n=541 over the timeseries). Clone E had thehighest selection coefficient (1 + s = 1.08 ± 0.043), having grown from undetectable proportions inearlier timepoints to 58% of cells by the end of the timeseries. Clone E also had the highest numberof breakpoints with 12.8 additional copy number breakpoints per cell, relative to the reference cloneC with the lowest (Figure A.10-c).71Figure 4.4: Comparison of fitness landscapes of breast cancer PDX models. (a-c) Clonal dynamics of HER2+ and TNBC models. (d)Heatmap representation of copy number profiles of 2,193 cells, grouped in 4 phylogenetic clades. (e) Phylogeny as per Figure 4.2 forHER2+. (f) Inferred fitClone trajectories and (g) selection coefficients for the HER2+ model. See Figure A.5 for the credible intervalsand Figure A.6 for pairwise comparison of the selection coefficients. (h-k) Analogous plots for the TNBC model (n=3,216 cells).724.3.2 Mixture experimentsIn the two previous sections, we used our framework to explore the evolutionary dynamics in threerelated cell lines and two independent PDX timeseries. Here we aim to examine the predictivecapabilities of our model. Specifically we are interested in verifying the selection coefficients es-timated by fitClone. As such we return to the following hypothesis: given a resolved startingconfiguration, and inferred selection coefficients in a sufficiently similar context, it is possible toestimate subclonal trajectories. Testing this hypothesis motivates an experimental design that wecall the mixture experiments. In this setup, we extract cells from an early (X3) and a late (X8)passage of the TNBC-SA609 series and then physically mix them. We established two such lines.In the first line (branch a) we aimed to have a mixture comprising equal proportions from the twotimepoints (a ration of 1.0 to 1.0). In the second line (branch b) we aimed to have 71.0% of the cellscome from the X3 passage and 29.0% of cells from the X8 passage (a ratio of 1.0 to 0.4). We did soto measure the effects of different starting population configurations on the selection coefficients.We call the initial physical mixture in each branch M0. For each branch, we describe the forecastingstudy and then the results from serially passaging the PDX line measured by DLP+.Mixture branch a We forward-simulated trajectories from fitClone using the median of theposterior distribution of the estimated selection coefficients (F = -0.01 ± 0.13, A = -0.00 ± 0.16, B= 0.00 ± 0.01, D = 0.00 ± 0.01, G = 0.02 ± 0.02, H = 0.03 ± 0.03, E = 0.08 ± 0.10) and starting-clonal proportions of (A=0.000 (no observed cells), B=0.07, C=0.25, D=0.51, E=0.02, F=0.000(no observed cells), G=0.08 H=0.07) Figures A.11-a, and 4.5-a,d. The starting clonal proportionswere inferred by adding the cells from the initial passage (M0) to the tree inferred from the cellsin the original series (TNBC-SA609) and assigning them to their respective clones. We generated10,000 trajectories from the model. Figure 4.5-c (top) shows the simulated trajectories in black.The mean clonal fraction at each step is shown in red. All clones except for clones E, D, and F arepredicted to vanish to clonal fractions of below 1%. We have combined the trajectories for clonesE and F since clone F (i) had fewer than 19 total cells in the original series and consequently itsselection coefficient had a high variance, (ii) is phylogenetically proximal to clone E (Figure 4.5-b)and thus likely represented a biologically similar population and (iii) finally, is not observed abovea threshold of 20 cells in any other line in the TNBC-SA609 family.73We experimentally tested these predictions by initiating a new PDX line with the remixedpopulation (from M0), serially passaged over 4 timepoints (Figure 4.5-a (top), and sequencedwith DLP+ (7,839 single cell genomes, median 1,354.5 per library)). After placing the cells fromall timepoints of this mixture experiment on the tree, they got assigned to seven clones fromthe original timeseries (all but clone A) with between 26 to 499 (median 162) cells (Figure 4.5-c(middle)). In Figure 4.5-c (top) blue dots show the observed clonal fractions at each timepointin PDX branch a. Clones with higher selection coefficients swept through the mixture timeseriesby passage 4. In the last timepoint, the clade composed of clones E and F comprised 94% ofcells, outcompeting low-fitness clones. The estimated selection coefficients were relatively stronglycorrelated (Pearson correlation of 0.795, considering only clones that reached overall prevalence ofover 1% in the original series, i.e., all clones except A and F).Mixture branch b We forward simulated 10,000 trajectories using identical selection coeffi-cient values as branch a, but different clonal proportions, estimated from adding cells from timepointX1 in mixture branch b to the TNBC-SA609 phylogenetic tree and assigning them to the corre-sponding clones (C = 0.02, D = 0.00, E = 0.05, F = 0.00, G = 0.06, H = 0.87). Unlike mixturebranch a, DLP+ data was not available for timepoint M0 in mixture branch b. We initiated aPDX line passaged over 5 timepoints that were sequenced using the DLP+ platform to generate6,730 single cell genomes (median 1,270 per library). We observed 6 clones from the original series,namely all clones except A and B. However, we only captured over the 5 timepoints, 1 cell fromclone D, 8 cells from clone F, and 13 cells from clone C. As predicted, clone E, which had thehighest predicted selection coefficient in the original series despite having started from a very lowclonal fraction (0.05), rises to an abundance of 0.28, while clone H steadily falls from 0.87 to about0.67.74allllllllllllllllllllllllllllllllllllllllllllllll lllllllllllllllllllll lllllllllllllllllllllllllllllllllllllllllllll lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll llllllB (26)C (84) D (162)E (499)F (268)G (148)H (496)Number of cellsl l l l l l l l0 1 7 20 43 79 130 198llllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll llllllllllllllllllllllllllllll llllllllllllllllllllllllllllllllllll lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll llllllll llllllllllllllllllllC (13)D (1)E (442)F (8)G (171)H (1891)Number of cellsl l l l l l l l0 1 7 20 43 79 130 198bll lllll ll l llllllll l ll ll llHD E_F GA B C0.00.10.20.00.10.20.00.10.20.00.40.80.00.40.80.00.40.8Diffusion time (Days/Ne)Clonal fraction0.000.250.500.751.000.00 0.05 0.10 0.15 0.20Diffusion time (days/Ne)Clonal fractionBDGHEF0.6 0.8 1.0 1.2 1.41 + s (relative to clone C)cll l l l ll lll ll lll ll l ll l l l lG HD E FA B C0.000.030.080.000.030.080.000.030.080.00.40.80.00.40.80.00.40.8Diffusion time (Days/Ne)Clonal fraction0.000.250.500.751.000.00 0.05 0.10 0.15 0.20Diffusion time (days/Ne)Clonal fractionFDHGE0.6 0.8 1.0 1.2 1.41 + s (relative to clone C)dFigure 4.5: Mixture experiments on TNBC-SA609 PDX. In each panel, mixture a is followed by mixture b. (a) Clonal proportions ofX3 and X8 to generate the initial mixture M0 and subsequent serial passaging, yielding 4 samples M1-M4 for mixture a and 5 samplesM1-M5 for mixture b. (b) Phylogenies showing cells observed in the mixture a (left) and mixture b (right) timeseries. (c) For mixturea: (left) Forward simulations using inferred selection coefficients and starting population proportions in the initial experimental mixture.Simulated trajectories are shown. (middle) Inferred trajectories of mixture timeseries. (right) Selection coefficients of fitClone fit toM1-M4 clonal abundance observations. See Figure A.5 for the credible intervals and Figure A.6 for pairwise comparison of the selectioncoefficients. (d) as in (c) but for mixture b. Note that in the prediction trajectories plot (Figure 4.5-c,d-(left)), the time-axis for theobservations is shrunk to best match the mean predicted trajectories line (red line).75Clone-specific fitness estimates forecast clonal competition trajectoriesThe analysis of the two mixture series presented here suggests that (i) we can validate the selectioncoefficients estimated from a timeseries using fitClone and (ii) it is possible make quantitativepredictions about the likely trajectories of tumour subpopulations at least at the clade level. Wenote that in the prediction trajectories plot (Figure 4.5-c,d-left), the time-axis for the observationsis shrunk to best match the mean predicted trajectories line (red line). The diffusion time horizonthat is obtained by dividing the generation time measured in days by the Ne estimate results intrajectories that are ahead of the biological system. While in branch a the original time horizon is0.24 and the best matching time is 0.20, in branch b the values are 0.23 and 0.08 respectively.4.4 In vivo systems under cisplatin treatmentSo far we have applied our framework to real-world datasets to track human cancer clones asidentified by their CN genotypes over time and to reason about their likely abundances in certainexperimental conditions. Now we investigate whether establishing baseline fitness measures couldhelp to interpret selection under drug administration. We first present the TNBC-SA609 series indetail and make an observation about early response to cisplatin treatment. We then analyse twoadditional TNBC PDX lines to check the reproduabiblty of our observation.4.4.1 TNBC-SA609 PDXWe tested how pharmacologic perturbation with cisplatin impacted the stability of the fitnesslandscape of the TNBC series. Using frozen material from the third timepoint (X3) in the originalTNBC-SA609 untreated line introduced in Section 4.3, we established a new line that was passaged4 times (X4, X5, X6, and X7) comprising 6,323 cells, 841 of which passed the QC step (Figure A.9).Figure A.1 shows how the different timepoints are related to one another. We also formed a separatebranch from the same starting material (frozen X3) and administered cisplatin (2mg/kg, Q3Dx8 i.p.max) serially over four successive passages to induce drug resistance (Figure A.11-b,c). For eachserially treated tumour in mice, a parallel set of transplanted mice were left untreated, establishingcorresponding drug ‘holiday’ samples (Figure 4.1-a). We coded the treated passages with ‘T’ anduntreated with ‘U’, initialised by the X3 untreated (U) passage. The first treatment passage (X476UT ) exhibited rapid tumour shrinkage (>50% of initial size). However X5 UTT, X6 UTTT andX7 UTTTT had progressively less response, indicating drug resistance and positive growth kinetics(Figure A.11-e). Decomposing the growth dynamics over (X3 U ; X4 UT ; X5 UTT ; X6 UTTT ; X7UTTTT ) into clonal trajectories with DLP+ analysis suggested that sustained cisplatin treatmentinverted the fitness landscape. A new clone R, derived from clone A in the phylogeny, but with adistinct clonal genotype (fewer copies of MYC and deletions at RB1, PRDM9 and NUDT15 loci(Figure 4.6-a,b)), swept to fixation comprising 48% (X4 UT ), 98% (X5 UTT ), 100% (X6 UTTT )and 100% (X7 UTTTT ) of cells across the treated series (Figure 4.6-c). Notably, the high fitnessclones E, H, G, D from the untreated series exhibited low fitness coefficients in the treatment seriesand were no longer detected (Figures 4.6-d, 4.9). Conversely, clones A, B, and C, comprising a lowfitness phylogenetic superclade, distinct from high fitness clones E and F in the untreated series,were the precursors to the resistant clone R (Figure 4.6-e). Thus, cisplatin perturbation resultedin a near complete reversal in the fitness landscape.Drug holidayNext, we asked whether the clonal dynamics in the presence of cisplatin were reversible by examiningthe drug holiday samples (Figure 4.1-a, X5 UTU ; X6 UTTU ; X7 UTTTU ). In the first drugholiday X5 UTU, clonal composition reverted to consist predominantly of precursor clone B with90% abundance, and only 10% abundance from clone R (Figure 4.9-(left)). However, in X6 UTTUand X7 UTTTU no reversion was detected, and these populations consisted of >99% clone R,similar to their on-treatment analogues. Thus, clonal competition in the absence of drug led toclones derived from the A-B clade outcompeting clone R, and clone-specific cisplatin resistancehas a fitness cost. Moreover, the genotype specificity of reversion between X4 UT to X5 UTUindicates that the clonal dynamics can be attributed to selection of genomically defined clones withdifferential fitness.Together, these data suggest the impact of cisplatin selective pressure on the starting tumourcell population is reversible while genomic clonal competition with precursor clones is still possible,but dominates the population once the evolutionary bottleneck narrows and purifies the population.77X5 UTTX5 UTUX7 UTTTTX7 UTTTUX6 UTTUX6 UTTTX3 UX4 UT+Rx-Rx+Rx-Rx+Rx-Rx+Rx-Rxcabd ePRDM9MYCRB1NUDT15Root+Rx-Rx+RxRootERRoot RootRootDCBBRBRR RRRRCBFigure 4.6: Impact of pharmacologic perturbation with cisplatin on fitness landscapes TNBC-SA609 PDX. (a) Copy number genotypeof clone E from untreated timeseries. (b) Copy number genotype of clone R from treated timeseries (arrows indicate differences to cloneE). (c) Clades with higher fitness in the untreated (-Rx) and treated (+Rx) series. (d) Evolution as a function of drug treatment anddrug holiday. Arrows indicate the path of serial passaging. For each sample, the phylogeny with clonal abundance from DLP+ is shown,reflecting selection. (e) Reversal in the fitness landscape - note clone A maps to clone R in post treatment since R is derived from A,but only emerges in the treatment series.784.4.2 Fitness landscape reversal in early response to cisplatinWe then asked if the reversal in the fitness landscape observed in the TNBC-SA609 system undercisplatin treatment is reproducible. To address this, we established two independent TNBC PDXsystems as follows.TNBC-SA535 We generated a total of 15,302 single cells out of which 4,023 passed our qualityfilters. TNBC-SA535 is a BRCA1 deficient patient derived tumour. We established a timeserieswith treated and untreated branches similar to TNBC-SA609 (see Figure A.2). We generated5 consecutively transplanted timepoints (X5, X6, X7, X8, X9) left untreated for a total of 1,341single cells (mean = 335, σ = 84.4 per timepoint). Simultaneously, we established a cisplatin treatedtimeseries starting from timepoint X6, and continued cisplatin treatment for 5 cycles up to X10,generating a total of 1,425 cells from scWGS (mean = 356, σ = 159 per timepoint) from 4 cycles. Acut of the phylogenetic tree inferred over all cells in this series, resulted in 7 clones. In the untreatedline, clonal fractions were A(0.003), B(0.702), C(0.034), D(0.006), E(0.006), F(0.013), and G(0.237).Clone B was chosen as the reference clone as it had a monotonically decreasing clonal fractiontrajectory in the untreated branch. Clonal trajectories were consistent with selection coefficientswith small relative differences in fitness (Figure 4.7-d). Clone G had the highest fitness (1 + s =1.01 ± 0.00751) closely followed by clone C (1 + s = 1.00 ± 0.0282) and the reference clone. In thetreated branch, clonal fractions were A(0.156), B(0.066), C(0.140), D(0.194), E(0.182), F(0.151),and G(0.112). In this regime, clone A emerged with the highest selection coefficient (1 + s = 1.03± 0.0152) followed by clone D (1 + s = 1.02 ± 0.0116). Notably clones A and D had low fitnessvalues under no treatment, whereas clones G (1+s = 1.01 ± 0.0115) and C (1+s = 1.02 ± 0.0119)had low fitness coefficients under treatment.TNBC-SA1035 Another independent PDX system with 14,170 single cells were generatedwhere 4,444 passed the quality filters. The experimental design diagram is shown in Figure A.3.An untreated branch with five serial passages (X4, X5, X6, X7, and X8) was established with atotal of 2,015 single cells. A parallel branch was treated with cisplatin starting at X5, X6, X7, andX8 comprising 1,596 filtered cells. 833 cells belonged to the drug-holiday timepoints and are notanalyzed here. Phylogenetic inference followed by cutting the tree yielded 11 clones (Figure 4.8).Clonal fractions over all timepoints in the untreated branch were A(0.097), B(0.140), C(0.087),79D(0.160), E(0.266), F(0.010), G(0.058), H(0.053), I(0.065), J(0.018), and K(0.047). Clone A’sabundance fell over time and it was chosen as the reference clone. Clone E rose from a clonal fractionof 0.028 at X4 to 0.69 at X8 and had the highest selection coefficient (1+s = 1.06 ± 0.0367). In thetreated branch, clonal fractions were A(0.065), B(0.129), C(0.132), D(0.066), E(0.055), F(0.018),G(0.205), H(0.144), I(0.094), J(0.014), and K(0.078). In this regime, G (1+s = 1.01 ± 0.0123) andH (1 + s = 1.02 ± 0.0135), which were among the clones with lower fitness in absence of treatment,rose to occupy 73% at X8 while clone E (1 + s = 0.993 ± 0.0344) fell from about 10% at X5 toundetectable at X8.80Figure 4.7: TNBC-SA535 PDX clonal dynamics with and without treatment. (a) Heatmap representation of copy number profiles of1,341 cells, grouped in 7 phylogenetic clades. (b) Phylogeny as per Figure 4.2 for TNBC-SA535 PDX untreated branch. (c) Observedclonal abundances. (d) Distribution over magnitude of difference between selection coefficients of pairs of clones. (e-h) Analogous plotsfor the treated branch (n=1,425 cells).81Figure 4.8: TNBC-SA1035 PDX clonal dynamics with and without treatment. (a) Heatmap representation of copy number profiles of2,015 cells, grouped in 11 phylogenetic clades. (b) Phylogeny as per Figure 4.2 for TNBC-SA1035 PDX untreated branch. (c) Observedclonal abundances. (d) Distribution over magnitude of difference between selection coefficients of pairs of clones. (e-h) Analogous plotsfor the treated branch (n=1,596 cells).82Figure 4.9 summarises the reversal in the fitness landscape in response to cisplatin treatment inTNBC-PDX model systems. In the TNBC-SA609 system the fitness landscape is inverted whereinclones more fit in the untreated regime (H, D) are less fit in the treated regime, whereas less fitclones in the untreated regime (A, B) are the most fit clones under treatment. This pattern ismirrored in two independent TNBC PDX lines treated with cisplatin, namely TNBC-SA535 andTNBC-SA1035. In TNBC-SA535, clones G, C, and B are drug-sensitive, while A and D are drug-resistant. The former have higher relative fitness in untreated versus treated regimes, while thelatter exhibit an inverted fitness pattern. Similarly, in TNBC-SA1035, drug-sensitive clones consistof clones E, B, and A, while the drug-resistant group comprises clones H, I, and G. From untreatedto treated, the first group goes from high to low fitness, while the second group goes from low tohigh fitness.83llllllllllllHBDGCAABCHGD−Rx +RxSA6090.000.250.500.751.000.0 0.2 0.4 0.6Diffusion time (days/Ne)Clonal fractionAGDBH0.6 0.8 1.0 1.2 1.41 + s (relative to clone C)Clones0.000.250.500.751.000.0 0.1 0.2 0.3 0.4Diffusion time (days/Ne)Clonal fractionDGHEFBA0.6 0.8 1.0 1.2 1.41 + s (relative to clone C)ClonesllllllllllllllGCBFDEAADFCEGB−Rx +RxSA5350.000.250.500.751.000.0 0.2 0.4 0.6Diffusion time (days/Ne)Clonal fractionAEDFCG0.6 0.8 1.0 1.2 1.41 + s (relative to clone B)Clones0.000.250.500.751.000.0 0.2 0.4 0.6Diffusion time (days/Ne)Clonal fractionGECFDA0.6 0.8 1.0 1.2 1.41 + s (relative to clone B)ClonesllllllllllllllllllllllEBADCHKJIGFHGICKBAEDFJ−Rx +RxSA10350.000.250.500.751.000.0 0.1 0.2 0.3 0.4Diffusion time (days/Ne)Clonal fractionFGIJKHCDBE0.6 0.8 1.0 1.2 1.41 + s (relative to clone A)Clones0.000.250.500.751.000.0 0.1 0.2 0.3 0.4Diffusion time (days/Ne)Clonal fractionJFDEBKCIGH0.6 0.8 1.0 1.2 1.41 + s (relative to clone A)ClonesFigure 4.9: Fitness landscape reversal in early cisplatin treatment in TNBC PDX models. We observe that clone specific resistance tocisplatin treatment arises in 3 independent TNBC PDX lines. In all three cases, clones with low fitness under no treatment exhibit highfitness under the treatment regime. In each panel, the left and right sub-panels are from the untreated and treated branches respectively.(top) Clones sorted by their median selection coefficient in -Rx and +Rx regimes. (middle) Inferred trajectories, and (bottom) selectioncoefficients of fitClone model fits to each branch. See Figure A.5 for the credible intervals and Figure A.6 for pairwise comparison ofthe selection coefficients.844.5 DiscussionHere we show that decoding the contributions of clonal competition landscapes in the course of tu-mour growth with and without drug selection can be achieved by measuring and modelling cellulardynamics with granular timeseries over several months-years. Single cell whole genome sequenc-ing allowed for the robust application of population genetic statistical models. TP53 ablation indiploid mammary epithelial cells resulted in clonal expansions driven by copy number alterations,with whole genome, whole chromosome and more localised alterations, all leading to fitness ad-vantages over diploid baseline. Positive selection attributed to copy number changes of all scalesmay be under-investigated [6]. This has implications for interpreting etiologic processes of TP53driven cancers where the rates of structural variation acquisition and deviation away from diploidconfigurations confer quantitative fitness advantages. In tumour evolution within patients, copynumber alterations, as a key biological process, is revealing impact on treatment outcomes andco-morbidities [105], innate and adaptive immune response [106, 107], the root cause of major ge-nomic reconfigurations [108] and evolutionary plasticity [16]. Our work here suggests that TP53mutation can lead to structural alterations in vitro that are also observed in breast cancer patients,thereby representing a realistic model for studying how the impact of driver mutations inducinggenomic instability leads to clonal expansions and evolutionary selection. Variously through singlecell approaches [109] or computational reconstruction of evolutionary histories [9], the evolutionaryimpact of structural variations is still a work in progress. Over successive generations in vitroand in PDX, in the context of pharmacologic and mutational perturbation, emergent copy numberchanges contribute to the kinetics of the fitness landscape, consistent with a continual diversifyingmechanism that induces competitive clonal advantages.The impact of drug intervention on cancer evolution is a key determinant of patient outcomesacross all human cancers. As clonal drug resistance was consistent with a fitness cost, we sug-gest this could be exploited in future therapeutic strategies. Forecasting the trajectories of cancerclones is of immediate importance to understanding therapeutic response in cancer and for de-ploying adaptive approaches [5]. The presence within a tumour of lineage precursors to resistantgenotypes may define time windows within which clonal competition could mediate plasticity totreatment. We suggest that population genetic modelling of timeseries tumour or tumour cell-free85DNA measurements to predict clonal evolution is tractable, but will require complementary mea-surements of genotypic clonal abundance to gain comprehensive understanding. Further study withtimeseries modelling will provide insight into therapeutic strategies promoting early intervention,drug combinations and evolution-aware approaches to clinical management [110].4.6 ConclusionsIn this chapter, we exemplified the computational methods that we developed in Chapters 2 and 3 oncell lines and PDX systems generated over a multi-year period. We used sitka and Lumberjack toidentify the clonal subpopulations in heterogenous single cell subpopulations and their abundancesover time based on their copy number profile. We then used fitClone to ascribe quantitativeselection coefficients to individual cancer clones and forecast competitive clonal dynamics over time.Specifically we used a mixture competition design to experimentally verify the selection coefficientsinferred by our model. We observed that fitness landscape is reversed in early response to treatmentwith cisplatin in three independent TNBC PDX timeseries. We hope that our results would helpthe research community map the fitness landscape of tumour development and ultimately aid inselecting better treatment regimes for patients with cancer.4.6.1 Limitations and future researchIn the context of pharmacological intervention, we are ultimately interested in mechanisms un-derlying sensitivity and resistance to drug. We acknowledge that establishing causal relationshipsis very difficult and requires controlled experiments and larger sample sizes. However, exploringthe potential correlates of differential fitness may aid in designing follow up experiments, and ineventually finding possible druggable molecular targets. In the TNBC-SA609 PDX series, we didnot find a single nucleotide variant (SNV) that was private to a clone and that was implicated inthe literature as oncogenic. It is possible that such a mutation exists but was missed either at se-quencing or by our computational pipelines. Follow up replication experiments sequenced at higherdepth may be required to rule out the existence of driving mutations. Other biological mechanismsincluding epigenetic changes may influence the fitness of clones and interrogating the system viaorthogonal data types such as chromatin accessibility and DNA methylation assays is the subject86of future research.In asexually reproducing populations, multiple subpopulations each with a set of beneficialmutations can coexist [111, 112]. TNBC-SA609 PDX series shows evidence of persistence of multipleclones over time 4.4-f. This competition between clones, called clonal interference, reduces the rateof adaptation and may lead to an underestimation of the selection coefficients of individual clones.In other words, in absence of competition from other clones, each clone may have a higher innateselection coefficient. Clonal interference can also result in more stochasticity in the evolutionaryprocess [112]. As a consequence, the clones detected in one timeseries, may behave differently inabsence of one or more of the original clones.Exponential growth is a good fit to the initial stage of tumour growth. Once the tumourhas grown to a certain size, it is expected that the growth would be boundary-driven. In anexponentially growing scenario, the observed clonal fraction trajectories more accurately representthe innate selection coefficients of clones, whereas in a boundary-driven growth scenario, spatiallimitations bias the inferred selection coefficients; for instance, a highly fit clone can emerge in themiddle of the tumour and due to space constraints will not grow as much a neutral clone growingat the boundary of the tumour [113].We note that in the formulation of fitClone, no pair of clones can be assigned an equal selectioncoefficient. As a consequence, fitClone cannot delineate clones that are identical from those thatare close in terms of evolutionary fitness. One remedy is to use a spike and slab prior that puts apoint mass on si = 0.87Chapter 5Conclusions5.1 Summary of contributionsIn this work we have presented a novel approach to investigate clonal evolution in longitudinalstudies in vitro and in vivo model systems. The methodological contributions include (i) the pre-processing, benchmarking and development of visualisations for a novel Bayesian phylogeneticsmethod (sitka) and its application in sorting single cells into clones in copy number space, alongwith an algorithm to cut the tree (Lumberjack), (ii) the adaptation of the Wright-Fisher diffusionmodel to track the trajectories of clonal fractions over time and (iii) the design and implementationof an inference engine to estimate evolutionary parameters and forecast in this model. We haveapplied these methods to investigate the predictability of clonal trajectories and quantify the re-sponse to cisplatin treatment in PDX models. Table 5.1 lists the repositories that host the sourcecode for the software implemented in this work.Name Description URL1 sitka-viz Visualisation suite forBayesian phylogeneticsvia sitka [1]https://github.com/UBC-Stat-ML/sitka-viz2 fitClone [2] The inference engine for theWright-Fisher model SDEhttps://github.com/UBC-Stat-ML/fitclone3 Lumberjack The tree-cutting algorithm forsitkahttps://github.com/UBC-Stat-ML/tree_cutting4 sitka-material Tools for preprocessing, post-processing, and benchmarkingof sitkahttps://github.com/UBC-Stat-ML/sitka-materialTable 5.1: List of repositories hosting the source code for software developed as part of this thesis.An important open question in the treatment of patients with cancer is why and how resistanceto therapy is developed. Advances in next generation single cell molecular sequencing as well aslongitudinal xenoengraftment has enabled us to more accurately interrogate and make quantitative88measurements of tumours. By serially transplanting a tumour sample into highly immunodeficientmice, a system called patient derived xenograft, one can continuously monitor how tumour compo-sition may evolve in a patient over time. This process generates timestamped samples comprisingthousands of cells that through SCS can be measured at thousands of genomic features includingCNA profiles. This is in contrast to the standard of care where often only up to two timepoints areacquired and measured at an aggregate course level. The new higher resolution datatype introducesmajor challenges in extracting meaningful knowledge from the data, including (i) how to identifybiologically meaningful groups of cells (i.e., clones) across multiple timepoints, and (ii) how to quan-titatively reason about the underlying evolutionary forces acting on the clones via their observeddynamics. We have developed tools to generate meaningful summaries, and analyse and visualisethousands of cells over thousands of CNAs. These include combining timeseries phylogenetic treesand copy number state plots, identifying clones, and inferring and comparing their evolutionaryfitness over time.To ascertain the evolutionary relationships of the cells comprising the tumour samples, a newBayesian phylogenetic inference method called sitka was developed in our research group. Wedeveloped a pipeline that accepts single cell CNAs as input, performs quality control and filteringdefined by the researcher, and transforms it into a binary space and a format ready for phylogeneticanalysis. The outputs of the analysis include the inferred consensus evolutionary tree, the predictedposterior binary genotype, and the error-corrected input matrices. We devised and implementedLumberjack, an algorithm that cuts a phylogenetic tree by grouping subsets of cells that aregenomically distinct based on both the topology of the tree and the cell’s CNAs. These subsets,called clones, are biologically important entities that are hypothesised to encode unique phenotypesincluding differential response to therapeutic interventions. The phylogenetic inference outputs arethen visualised in two modes, detailed and publication quality. In the former, the tree is alignedand plotted next to a heatmap representation of the CNAs displaying the genome of individualcells and their placement on the tree, interlaced with the posterior binary genotype state and theerror-corrected input matrices. This has helped in interpreting and inspecting the inference results.The latter mode produces summarised trees organised in time steps reflective of the longitudinalexperimental designs where changes in the tumour composition can be immediately gleaned via thechange in colour (clone assignment) and size (number of constituent cells) of the nodes on the tree.89We used this visualisation to help characterise a unique mode of acquired reversible drug-resistancein breast cancers in response to cisplatin resistance.Tracking the relative abundances of the clones over time produces a timeseries reminiscent ofallele frequencies. We adopted ideas from the population genetics literature to model the evolutionof tumour growth via a state space model, the dynamics of which is dictated by the diffusion ap-proximation to the Wright-Fisher process. This framework, called fitClone induces an intractablelikelihood that makes inference difficult. To address this, we extended a state of the art pseudo-marginal sequential Monte Carlo algorithm and implemented the inference engine. Briefly, thealgorithm infers evolutionary fitness coefficients for each clone by simulating trajectories that areconsistent with their observed clonal abundances. We devised a novel algorithmic extension thatscaled simultaneous inference to dozens of clones. We then detected computational bottlenecks andre-implemented them in an efficient vectorised manner.5.2 Future research directionsThe Wright-Fisher model investigated in this work assumes all clones exist at the initial timepoint.If their observed clonal fraction is zero at the first timepoint, the model assumes that the clone existsat fractions under the threshold of detection. This ignores the possibility of de novo mutations andclones that arise only at later timepoints. One drawback is that estimated selection coefficientswould be biased as the model has to account for the absence and then sudden rise of such clonesby decreasing the selection coefficient. One potential remedy is to incorporate mutations using theWright-Fisher model with selection and mutation.The selection coefficient of each clone is also assumed to be constant over time. A step modelthat allows for change in the selection coefficient may produce a better fit to the data. The sitkaphylogenetics model does not infer branch lengths. One potential direction is to use the time-stamps of the single cells under each edge to infer a relative evolutionary change rate, that is tocreate a clocked model.In this work we are restricted to measurements of clonal fractions over time, and are agnosticto the tumour microenvironment (TME), the spatial placement of clones, their epigenetic profile,protein levels, and clinical covariates, to name a few. Cancerous cells are situated in their TME and90their interactions can shape tumour fate through mechanisms such as immune escape [114–116].More broadly other evolutionary paradigms, clonal competition, cooperation, etc, are not investi-gated here. Furthermore we do not take into account any medical correlates including patient’sage, smoking status, and family history. With the availability of larger longitudinal datasets, itmay be possible to elucidate the effects of these covariates, perhaps via embedding the fitClonemodel in the framework of generalised linear models.In Chapter 4 we identified clones using their genomic CN profiles. Establishing what genotypicfeatures incur higher fitness is necessary for reliable rational treatment strategies. Epistasis in whichmultiple mutations across multiple genomic loci interact to control a phenotype are more difficultto identify. Further controlled experiments in much larger cohorts could be used to characterisesuch relationships.In this work we only considered biological substrates acquired from solid tumours. Difficultiesof obtaining biopsies from patients may hinder continual monitoring of tumour growth and clonalcomposition of patients over time. One potential way to alleviate this is to test for circulatingtumour DNA (ctDNA) in blood samples. It has been shown that ctDNA has good fidelity to theoriginal tumour in some cancers [117, 118]. It would be interesting to test the robustness andapplicability of our framework to this type of data.Advances in non-invasive imaging [119] establish a localised view at micro (histopathology)and macro (radiology) levels. Most diagnosis and prognosis routines involve the use of one familyof data, e.g., one or more molecular-assays including genomics/transcriptomics/proteomics. Someeffort has been made to leverage more than one data family, often combining imaging and molecularassays [120–122]. This line of research has established the complementary nature of multi-modaldata, where for instance multiple genomic states of tumours such as aneuploidy are highly correlatedwith quantitative features extracted from histopathological images, and resulted in characterisingsome aspects of the tumour-microenvironment interaction [114]. While these models outperformdiagnosis solely based on one family of datasets, the gain over the standard of care, typicallyinvolving grading and staging by a pathologist based on histology slides and input from molecularassays is negligible [121]. An interesting research direction is exploring the integrative causal multi-modal frameworks which may improve the accuracy of diagnosis and prognosis [123–126].915.3 Concluding remarksCancer is a complex phenomenon where the laws of Darwinian evolution can lead to a uniquedisease in each patient. This calls for a personalised rational treatment regime where the spe-cific genomic aberrations of the patient’s cancer are taken into account. Advancements in scalablemolecular sequencing technologies allow for more accurate and more frequent monitoring of tumourdevelopment. Also, in the case of solid tumours where obtaining biopsies could be difficult, usingctDNA instead shows promise as a less invasive method to continuously track the clonal makeupof tumours. It can then be determined whether a tumour is relapsing and if so, which constitutingclone(s) is responsible for the resurgence. Computational methods that model the dynamics oftumour development can be employed to prioritize the target and timing for treatment admin-istration. The body of work developed here has been used in a study [2] to elucidate the earlymechanism of tumour resistance to cisplatin. We hope that with availability of more longitudinaldata, our framework would help the research community to map the fitness landscape of tumourdevelopment with and without intervention. This could aid in selecting better treatment regimesfor patients with cancer.92Bibliography[1] Fatemeh Dorri, Sohrab Salehi, Kevin Chern, Tyler Funnell, Marc Williams, Daniel Lai,Mirela Andronescu, Kieran R Campbell, Andrew McPherson, Samuel Aparicio, Andrew Roth,Sohrab P Shah, and Alexandre Bouchard-Coˆte´. Efficient Bayesian inference of phylogenetictrees from large scale, low-depth genome-wide single-cell data. bioRxiv, 2020.[2] S. Salehi, F. Kabir, N. Ceglia, M. Andronescu, M. Williams, K. R. Campbell, T. Ma-sud, B. Wang, J. Biele, J. Brimhal, J. Ting, A. W. Zhang, C. O’Flanagan, F. Dorri,N. Rusk, E. Laks, H. Lee, T. Algara, S. Lee, B. Y. C. Cheng, P. Eirew, T. Kono, J. Pham,D. Grewel, D. Lai, R. Moore, A. J. Mungall, M. A Marra, IMAXT Consertium, A. McPher-son, A. Bouchard-Coˆte´, S. Aparicio, and S. P. Shah. Single cell fitness landscapes induced bygenetic and pharmacologic perturbations in cancer. To be published, 2020.[3] Sydney M Shaffer, Margaret C Dunagin, Stefan R Torborg, Eduardo A Torre, BenjaminEmert, Clemens Krepler, Marilda Beqiri, Katrin Sproesser, Patricia A Brafford, Min Xiao,Elliott Eggan, Ioannis N Anastopoulos, Cesar A Vargas-Garcia, Abhyudai Singh, Kather-ine L Nathanson, Meenhard Herlyn, and Arjun Raj. Rare cell variability and drug-inducedreprogramming as a mode of cancer drug resistance. Nature, 546(7658):431–435, June 2017.[4] Chong Sun, Liqin Wang, Sidong Huang, Guus Jje Heynen, Anirudh Prahallad, CarolineRobert, John Haanen, Christian Blank, Jelle Wesseling, Stefan M Willems, and Others.Reversible and adaptive resistance to BRAF (V600E) inhibition in melanoma. Nature,508(7494):118–122, 2014.[5] Neil Vasan, Jose´ Baselga, and David M Hyman. A view on drug resistance in cancer. Nature,575(7782):299–309, November 2019.[6] Devon A. Lukow, Erin L. Sausville, Pavit Suri, Narendra Kumar Chunduri, Justin Leu, Jude93Kendall, Zihua Wang, Zuzana Storchova, and Jason M. Sheltzer. Chromosomal instabilityaccelerates the evolution of resistance to anti-cancer therapies. bioRxiv, 2020.[7] Marc J Williams, Benjamin Werner, Chris P Barnes, Trevor A Graham, and Andrea Sottoriva.Identification of neutral tumor evolution across cancer types. Nature genetics, 48(3):238, 2016.[8] Marc J Williams, Benjamin Werner, Timon Heide, Christina Curtis, Chris P Barnes, AndreaSottoriva, and Trevor A Graham. Quantification of subclonal selection in cancer from bulksequencing data. Nature genetics, page 1, 2018.[9] Moritz Gerstung, Clemency Jolly, Ignaty Leshchiner, Stefan C Dentro, Santiago Gonza-lez, Daniel Rosebrock, Thomas J Mitchell, Yulia Rubanova, Pavana Anur, Kaixian Yu,Maxime Tarabichi, Amit Deshwar, Jeff Wintersinger, Kortine Kleinheinz, Ignacio Va´zquez-Garc´ıa, Kerstin Haase, Lara Jerman, Subhajit Sengupta, Geoff Macintyre, Salem Malikic,Nilgun Donmez, Dimitri G Livitz, Marek Cmero, Jonas Demeulemeester, Steven Schumacher,Yu Fan, Xiaotong Yao, Juhee Lee, Matthias Schlesner, Paul C Boutros, David D Bowtell,Hongtu Zhu, Gad Getz, Marcin Imielinski, Rameen Beroukhim, S Cenk Sahinalp, Yuan Ji,Martin Peifer, Florian Markowetz, Ville Mustonen, Ke Yuan, Wenyi Wang, Quaid D Morris,PCAWG Evolution & Heterogeneity Working Group, Paul T Spellman, David C Wedge, Pe-ter Van Loo, and PCAWG Consortium. The evolutionary history of 2,658 cancers. Nature,578(7793):122–128, February 2020.[10] S P Shah, A Roth, R Goya, A Oloumi, G Ha, Y Zhao, G Turashvili, J Ding, K Tse, G Haffari,A Bashashati, L M Prentice, J Khattra, A Burleigh, D Yap, V Bernard, A McPherson,K Shumansky, A Crisan, R Giuliany, A Heravi-Moussavi, J Rosner, D Lai, I Birol, R Varhol,A Tam, N Dhalla, T Zeng, K Ma, S K Chan, M Griffith, A Moradian, S W Cheng, G B Morin,P Watson, K Gelmon, S Chia, S F Chin, C Curtis, O M Rueda, P D Pharoah, S Damaraju,J Mackey, K Hoon, T Harkins, V Tadigotla, M Sigaroudinia, P Gascard, T Tlsty, J F Costello,I M Meyer, C J Eaves, W W Wasserman, S Jones, D Huntsman, M Hirst, C Caldas, M AMarra, and S Aparicio. The clonal and mutational evolution spectrum of primary triple-negative breast cancers. Nature, 486(7403):395–399, June 2012.[11] S Nik-Zainal, P Van Loo, D C Wedge, L B Alexandrov, C D Greenman, K W Lau, K Raine,94D Jones, J Marshall, M Ramakrishna, and Others. The life history of 21 breast cancers. Cell,149(5):994–1007, 2012.[12] In˜igo Martincorena, Keiran M Raine, Moritz Gerstung, Kevin J Dawson, Kerstin Haase,Peter Van Loo, Helen Davies, Michael R Stratton, and Peter J Campbell. Universal patternsof selection in cancer and somatic tissues, 2017.[13] Khurum H Khan, David Cunningham, Benjamin Werner, Georgios Vlachogiannis, Inmac-ulada Spiteri, Timon Heide, Javier Fernandez Mateos, Alexandra Vatsiou, Andrea Lampis,Mahnaz Darvish Damavandi, Hazel Lote, Ian Said Huntingford, Somaieh Hedayat, Ian Chau,Nina Tunariu, Giulia Mentrasti, Francesco Trevisani, Sheela Rao, Gayathri Anandappa,David Watkins, Naureen Starling, Janet Thomas, Clare Peckitt, Nasir Khan, Massimo Rugge,Ruwaida Begum, Blanka Hezelova, Annette Bryant, Thomas Jones, Paula Proszek, MatteoFassan, Jens C Hahne, Michael Hubank, Chiara Braconi, Andrea Sottoriva, and Nicola Va-leri. Longitudinal liquid biopsy and mathematical modeling of clonal evolution forecast timeto treatment failure in the PROSPECT-C phase II colorectal cancer clinical trial, 2018.[14] Marco Gerlinger, Stuart Horswell, James Larkin, Andrew J Rowan, Max P Salm, IgnacioVarela, Rosalie Fisher, Nicholas McGranahan, Nicholas Matthews, Claudio R Santos, PierreMartinez, Benjamin Phillimore, Sharmin Begum, Adam Rabinowitz, Bradley Spencer-Dene,Sakshi Gulati, Paul A Bates, Gordon Stamp, Lisa Pickering, Martin Gore, David L Nicol,Steven Hazell, P Andrew Futreal, Aengus Stewart, and Charles Swanton. Genomic architec-ture and evolution of clear cell renal cell carcinomas defined by multiregion sequencing. Nat.Genet., 46(3):225–233, March 2014.[15] Mariam Jamal-Hanjani, Gareth A Wilson, Nicholas McGranahan, Nicolai J Birkbak,Thomas BK Watkins, Selvaraju Veeriah, Seema Shafi, Diana H Johnson, Richard Mitter,Rachel Rosenthal, et al. Tracking the evolution of non–small-cell lung cancer. New EnglandJournal of Medicine, 376(22):2109–2121, 2017.[16] Saioa Lo´pez, TRACERx Consortium, Emilia L Lim, Stuart Horswell, Kerstin Haase, ArianaHuebner, Michelle Dietzen, Thanos P Mourikis, Thomas B K Watkins, Andrew Rowan,Sally M Dewhurst, Nicolai J Birkbak, Gareth A Wilson, Peter Van Loo, Mariam Jamal-95Hanjani, Charles Swanton, and Nicholas McGranahan. Interplay between whole-genomedoubling and the accumulation of deleterious alterations in cancer evolution, 2020.[17] Andrew McPherson, Andrew Roth, Emma Laks, Tehmina Masud, Ali Bashashati, Allen WZhang, Gavin Ha, Justina Biele, Damian Yap, Adrian Wan, Leah M Prentice, JaswinderKhattra, Maia A Smith, Cydney B Nielsen, Sarah C Mullaly, Steve Kalloger, AnthonyKarnezis, Karey Shumansky, Celia Siu, Jamie Rosner, Hector Li Chan, Julie Ho, NataliyaMelnyk, Janine Senz, Winnie Yang, Richard Moore, Andrew J Mungall, Marco A Marra,Alexandre Bouchard-Coˆte´, C Blake Gilks, David G Huntsman, Jessica N McAlpine, SamuelAparicio, and Sohrab P Shah. Divergent modes of clonal spread and intraperitoneal mixingin high-grade serous ovarian cancer. Nat. Genet., 48(7):758–767, 2016.[18] Caroline J Watson, A L Papula, Gladys Y P Poon, Wing H Wong, Andrew L Young, Todd EDruley, Daniel S Fisher, and Jamie R Blundell. The evolutionary dynamics and fitnesslandscape of clonal hematopoiesis, 2020.[19] Benjamin H Good, Michael J McDonald, Jeffrey E Barrick, Richard E Lenski, and Michael MDesai. The dynamics of molecular evolution over 60,000 generations. Nature, 551(7678):45–50,November 2017.[20] Emma Laks, Andrew McPherson, Hans Zahn, Daniel Lai, Adi Steif, Jazmine Brimhall,Justina Biele, Beixi Wang, Tehmina Masud, Jerome Ting, Diljot Grewal, Cydney Nielsen,Samantha Leung, Viktoria Bojilova, Maia Smith, Oleg Golovko, Steven Poon, Peter Eirew,Farhia Kabeer, Teresa Ruiz de Algara, So Ra Lee, M Jafar Taghiyar, Curtis Huebner, JessicaNgo, Tim Chan, Spencer Vatrt-Watts, Pascale Walters, Nafis Abrar, Sophia Chan, MattWiens, Lauren Martin, R Wilder Scott, T Michael Underhill, Elizabeth Chavez, Chris-tian Steidl, Daniel Da Costa, Yussanne Ma, Robin J N Coope, Richard Corbett, StephenPleasance, Richard Moore, Andrew J Mungall, Colin Mar, Fergus Cafferty, Karen Gel-mon, Stephen Chia, CRUK IMAXT Grand Challenge Team, Marco A Marra, Carl Hansen,Sohrab P Shah, and Samuel Aparicio. Clonal decomposition and DNA replication states de-fined by scaled Single-Cell genome sequencing. Cell, 179(5):1207–1221.e22, November 2019.[21] Hans Zahn, Adi Steif, Emma Laks, Peter Eirew, Michael VanInsberghe, Sohrab P Shah,96Samuel Aparicio, and Carl L Hansen. Scalable whole-genome single-cell library preparationwithout preamplification. Nature methods, 14(2):167, 2017.[22] Andrew Roth, Jaswinder Khattra, Damian Yap, Adrian Wan, Emma Laks, Justina Biele,Gavin Ha, Samuel Aparicio, Alexandre Bouchard-Coˆte´, and Sohrab P Shah. Pyclone: sta-tistical inference of clonal population structure in cancer. Nature methods, 11(4):396–398,2014.[23] Sohrab Salehi, Adi Steif, Andrew Roth, Samuel Aparicio, Alexandre Bouchard-Coˆte´, andSohrab P Shah. ddClone: joint statistical inference of clonal populations from single cell andbulk tumour sequencing data. Genome biology, 18(1):44, 2017.[24] Yi Kan Wang, Ali Bashashati, Michael S Anglesio, Dawn R Cochrane, Diljot S Grewal, GavinHa, Andrew McPherson, Hugo M Horlings, Janine Senz, Leah M Prentice, et al. Genomicconsequences of aberrant DNA repair mechanisms stratify ovarian cancer histotypes. Naturegenetics, 49(6):856, 2017.[25] Caitriona Holohan, Sandra Van Schaeybroeck, Daniel B Longley, and Patrick G Johnston.Cancer drug resistance: an evolving paradigm. Nature Reviews Cancer, 13(10):714, 2013.[26] Xueli Nan, Chao Xie, Xueyan Yu, and Jie Liu. EGFR TKI as first-line treatment for patientswith advanced EGFR mutation-positive non-small-cell lung cancer. Oncotarget, 8(43):75712,2017.[27] Sandra Misale, Rona Yaeger, Sebastijan Hobor, Elisa Scala, Manickam Janakiraman, DavidLiska, Emanuele Valtorta, Roberta Schiavo, Michela Buscarino, Giulia Siravegna, et al. Emer-gence of KRAS mutations and acquired resistance to anti-EGFR therapy in colorectal cancer.Nature, 486(7404):532, 2012.[28] William A Hammond, Abhisek Swaika, and Kabir Mody. Pharmacologic resistance in col-orectal cancer: a review. Therapeutic advances in medical oncology, 8(1):57–84, 2016.[29] Tomas Reinert, Everardo D Saad, Carlos H Barrios, and Jose´ Bines. Clinical implications ofESR1 mutations in hormone receptor-positive advanced breast cancer. Frontiers in oncology,7:26, 2017.97[30] Lindsay Angus, Nick Beije, Agnes Jager, John WM Martens, and Stefan Sleijfer. ESR1mutations: Moving towards guiding treatment decision-making in metastatic breast cancerpatients. Cancer treatment reviews, 52:33–40, 2017.[31] Benjamin G Bitler, Zachary L Watson, Lindsay J Wheeler, and Kian Behbakht. PARPinhibitors: clinical utility and possibilities of overcoming resistance. Gynecologic oncology,2017.[32] Dan Graur and Wen-Hsiung Li. Fundamentals of Molecular Evolution. Sinauer AssociatesInc., 2nd edition, 2000.[33] Haiwei Luo and Austin L Hughes. dN/dS does not show positive selection drives separationof polar-tropical SAR11 populations. Molecular systems biology, 8(1):625, 2012.[34] Malachi Griffith, Christopher A Miller, Obi L Griffith, Kilannin Krysiak, Zachary L Skidmore,Avinash Ramu, Jason R Walker, Ha X Dang, Lee Trani, David E Larson, et al. Optimizingcancer genome sequencing and analysis. Cell systems, 1(3):210–223, 2015.[35] Andrea Sottoriva, Haeyoun Kang, Zhicheng Ma, Trevor A Graham, Matthew P Salomon,Junsong Zhao, Paul Marjoram, Kimberly Siegmund, Michael F Press, Darryl Shibata, et al.A Big Bang model of human colorectal tumor growth. Nature genetics, 47(3):209, 2015.[36] Sarah P Otto and Troy Day. A biologist’s guide to mathematical modeling in ecology andevolution. Princeton University Press, 2011.[37] Jason Moffat, Dorre A Grueneberg, Xiaoping Yang, So Young Kim, Angela M Kloepfer,Gregory Hinkle, Bruno Piqani, Thomas M Eisenhaure, Biao Luo, Jennifer K Grenier, et al.A lentiviral RNAi library for human and mouse genes applied to an arrayed viral high-contentscreen. Cell, 124(6):1283–1298, 2006.[38] Ophir Shalem, Neville E Sanjana, Ella Hartenian, Xi Shi, David A Scott, Tarjei S Mikkelsen,Dirk Heckl, Benjamin L Ebert, David E Root, John G Doench, et al. Genome-scale CRISPR-Cas9 knockout screening in human cells. Science, 343(6166):84–87, 2014.98[39] Bernhard Schmierer, Sandeep K Botla, Jilin Zhang, Mikko Turunen, Teemu Kivioja, andJussi Taipale. CRISPR/Cas9 screening using unique molecular identifiers. Molecular systemsbiology, 13(10):945, 2017.[40] Zoe¨ N Rogers, Christopher D McFarland, Ian P Winters, Santiago Naranjo, Chen-HuaChuang, Dmitri Petrov, and Monte M Winslow. A quantitative and multiplexed approach touncover the fitness landscape of tumor suppression in vivo. Nature methods, 14(7):737, 2017.[41] Thomas O. McDonald, Shaon Chakrabarti, and Franziska Michor. Currently available bulksequencing data do not necessarily support a model of neutral tumor evolution. NatureGenetics, 2018.[42] Timon Heide, Luis Zapata, Marc J. Williams, Benjamin Werner, Giulio Caravagna, Chris P.Barnes, Trevor A. Graham, and Andrea Sottoriva. Reply to ‘Neutral tumor evolution?’.Nature Genetics, 2018.[43] Peter Eirew, Adi Steif, Jaswinder Khattra, Gavin Ha, Damian Yap, Hossein Farahani, KarenGelmon, Stephen Chia, Colin Mar, Adrian Wan, et al. Dynamics of genomic clones in breastcancer patient xenografts at single-cell resolution. Nature, 518(7539):422–426, 2015.[44] Juliet A Williams. Using PDX for preclinical cancer drug discovery: the evolving field.Journal of clinical medicine, 7(3):41, 2018.[45] Christopher D Willey, Ashley N Gilbert, Joshua C Anderson, and George Yancey Gillespie.Patient-derived xenografts as a model system for radiation research. In Seminars in radiationoncology, volume 25, pages 273–280. Elsevier, 2015.[46] Gaelen T Hess, Laure Fre´sard, Kyuho Han, Cameron H Lee, Amy Li, Karlene A Cimprich,Stephen B Montgomery, and Michael C Bassik. Directed evolution using dCas9-targetedsomatic hypermutation in mammalian cells. Nature methods, 13(12):1036, 2016.[47] Kaston Leung, Anders Klaus, Bill K Lin, Emma Laks, Justina Biele, Daniel Lai, AliBashashati, Yi-Fei Huang, Radhouane Aniba, Michelle Moksa, et al. Robust high-performance nanoliter-volume single-cell multiple displacement amplification on planar sub-strates. Proceedings of the National Academy of Sciences, 113(30):8484–8489, 2016.99[48] Yong Wang, Jill Waters, Marco L Leung, Anna Unruh, Whijae Roh, Xiuqing Shi, Ken Chen,Paul Scheet, Selina Vattathil, Han Liang, et al. Clonal evolution in breast cancer revealed bysingle nucleus genome sequencing. Nature, 512(7513):155, 2014.[49] E. M. Ross and F. Markowetz. OncoNEM: inferring tumor evolution from single-cell sequenc-ing data. Genome Biology, 17:69, 2016.[50] Russell Schwartz and Alejandro A Scha¨ffer. The evolution of tumour phylogenetics: principlesand practice. Nature Reviews Genetics, 18(4):213–229, 2017.[51] Katharina Jahn, Jack Kuipers, and Niko Beerenwinkel. Tree inference for single-cell data.Genome biology, 17(1):1–17, 2016.[52] Christopher A Miller, Brian S White, Nathan D Dees, Malachi Griffith, John S Welch, Obi LGriffith, Ravi Vij, Michael H Tomasson, Timothy A Graubert, Matthew J Walter, et al.SciClone: inferring clonal architecture and tracking the spatial and temporal patterns oftumor evolution. PLoS Comput Biol, 10(8):e1003665, 2014.[53] J. Singer, J. Kuipers, K. Jahn, and N. Beerenwinkel. SCIΦ: Single-cell mutation identificationvia phylogenetic inference. bioRxiv, page 290908, March 2018.[54] Fang Wang, Qihan Wang, Vakul Mohanty, Shaoheng Liang, Jinzhuang Dou, Jincheng Han,Darlan Conterno Minussi, Ruli Gao, Li Ding, Nicholas Navin, and Ken Chen. Single-cell copynumber lineage tracing enabling gene discovery. bioRxiv, 2020.[55] X. Xu, Y. Hou, X. Yin, L. Bao, A. Tang, L. Song, F. Li, S. Tsang, K. Wu, H. Wu, W. He,L. Zeng, M. Xing, R. Wu, H. Jiang, X. Liu, D. Cao, G. Guo, X. Hu, Y. Gui, Z. Li, W. Xie,X. Sun, M. Shi, Z. Cai, B. Wang, M. Zhong, J. Li, Z. Lu, N. Gu, X. Zhang, L. Goodman,L. Bolund, J. Wang, H. Yang, K. Kristiansen, M. Dean, Y. Li, and J. Wang. Single-CellExome Sequencing Reveals Single-Nucleotide Mutation Characteristics of a Kidney Tumor.Cell, 148(5):886–895, March 2012.[56] T. L. Williams and B. M. E. Moret. An investigation of phylogenetic likelihood methods.In Third IEEE Symposium on Bioinformatics and Bioengineering, 2003. Proceedings., pages79–86, March 2003.100[57] D.B. Wilson. Generating Random Spanning Trees More Quickly Than the Cover Time. InProceedings of the Twenty-eighth Annual ACM Symposium on Theory of Computing, STOC’96, pages 296–303, New York, NY, USA, 1996. ACM.[58] Radford M Neal. Slice sampling. Annals of statistics, pages 705–741, 2003.[59] A. Bouchard-Coˆte´, K. Chern, D. Cubranic, S. Hosseini, J. Hume, M. Lepur, Z. Ouyang,and G. Sgarbi. Blang: Bayesian declarative modelling of arbitrary data structures.arXiv:1912.10396 [stat], 2019.[60] Christian Robert. The Bayesian choice: from decision-theoretic foundations to computationalimplementation. Springer Science & Business Media, 2007.[61] William J Youden. Index for rating diagnostic tests. Cancer, 3(1):32–35, 1950.[62] Alan E Gelfand. Model determination using sampling-based methods. Markov chain MonteCarlo in practice, pages 145–161, 1996.[63] R.R. Sokal, C.D. Michener, and University of Kansas. A Statistical Method for EvaluatingSystematic Relationships. University of Kansas science bulletin. University of Kansas, 1958.[64] Naruya Saitou and Masatoshi Nei. The neighbor-joining method: a new method for recon-structing phylogenetic trees. Molecular biology and evolution, 4(4):406–425, 1987.[65] Leland McInnes, John Healy, and James Melville. Umap: Uniform manifold approximationand projection for dimension reduction. arXiv preprint arXiv:1802.03426, 2018.[66] Ricardo JGB Campello, Davoud Moulavi, and Jo¨rg Sander. Density-based clustering basedon hierarchical density estimates. In Pacific-Asia conference on knowledge discovery and datamining, pages 160–172. Springer, 2013.[67] John P Huelsenbeck and Fredrik Ronquist. MRBAYES: Bayesian inference of phylogenetictrees. Bioinformatics, 17(8):754–755, 2001.[68] Daniel L Ayres, Michael P Cummings, Guy Baele, Aaron E Darling, Paul O Lewis, David LSwofford, John P Huelsenbeck, Philippe Lemey, Andrew Rambaut, and Marc A Suchard.101Beagle 3: improved performance, scaling, and usability for a high-performance computinglibrary for statistical phylogenetics. Systematic biology, 68(6):1052–1061, 2019.[69] Kijong Yi and Young Seok Ju. Patterns and mechanisms of structural variations in humancancer. Experimental & Molecular Medicine, 50(8):98, 2018.[70] Sweta Mishra and Johnathan R. Whetstine. Different facets of copy number changes: Per-manent, transient, and adaptive. Molecular and Cellular Biology, 36(7):1050–1063, 2016.[71] Richard A Neher and Oskar Hallatschek. Genealogies of rapidly adapting populations. Pro-ceedings of the National Academy of Sciences, 110(2):437–442, 2013.[72] Jeffrey J. Quinn, Matthew G. Jones, Ross A. Okimoto, Shigeki Nanjo, Michelle M. Chan,Nir Yosef, Trever G. Bivona, and Jonathan S. Weissman. Single-cell lineages reveal the rates,routes, and drivers of metastasis in cancer xenografts. Science, 2021.[73] Jeffrey E Barrick, Dong Su Yu, Sung Ho Yoon, Haeyoung Jeong, Tae Kwang Oh, DominiqueSchneider, Richard E Lenski, and Jihyun F Kim. Genome evolution and adaptation in along-term experiment with Escherichia coli. Nature, 461(7268):1243–1247, 2009.[74] Jose Blanchet. Exact simulation for multivariate Itoˆ diffusions. 2017.[75] Paula Tataru, Maria Simonsen, Thomas Bataillon, and Asger Hobolth. Statistical inferencein the Wright–Fisher model using allele frequency data. Systematic biology, 66(1):e30–e46,2017.[76] Mark A Beaumont, Jean-Marie Cornuet, Jean-Michel Marin, and Christian P Robert. Adap-tive approximate Bayesian computation. Biometrika, 96(4):983–990, 2009.[77] Matthieu Foll, Yu-Ping Poh, Nicholas Renzette, Anna Ferrer-Admetlla, Claudia Bank, Hyun-jin Shim, Anna-Sapfo Malaspinas, Gregory Ewing, Ping Liu, Daniel Wegmann, et al. In-fluenza virus drug resistance: a time-sampled population genetics perspective. PLoS Genet,10(2):e1004185, 2014.[78] Jonathan P Bollback, Thomas L York, and Rasmus Nielsen. Estimation of 2Nes from temporalallele frequency data. Genetics, 179(1):497–502, 2008.102[79] Anna-Sapfo Malaspinas, Orestis Malaspinas, Steven N Evans, and Montgomery Slatkin. Es-timating allele age and selection coefficient from time-serial data. Genetics, 192(2):599–607,2012.[80] Anna Ferrer-Admetlla, Christoph Leuenberger, Jeffrey D Jensen, and Daniel Wegmann. Anapproximate markov model for the Wright–Fisher diffusion and its application to time seriesdata. Genetics, 203(2):831–846, 2016.[81] Alexandros Beskos, Gareth O Roberts, et al. Exact simulation of diffusions. The Annals ofApplied Probability, 15(4):2422–2444, 2005.[82] Murray Pollock, Adam M Johansen, Gareth O Roberts, et al. On the exact and ε-strongsimulation of (jump) diffusions. Bernoulli, 22(2):794–856, 2016.[83] Paul A Jenkins, Dario Spano, et al. Exact simulation of the Wright–Fisher diffusion. TheAnnals of Applied Probability, 27(3):1478–1509, 2017.[84] Simon JA Malham and Anke Wiese. An introduction to SDE simulation. arXiv preprintarXiv:1004.0646, 2010.[85] Fredrik Lindsten, Michael I Jordan, and Thomas B Scho¨n. Particle Gibbs with ancestorsampling. The Journal of Machine Learning Research, 15(1):2145–2184, 2014.[86] Fredrik Lindsten, Pete Bunch, Sumeetpal S Singh, and Thomas B Scho¨n. Particle an-cestor sampling for near-degenerate or intractable state transition models. arXiv preprintarXiv:1505.06356, 2015.[87] PE Jorde, S Palm, and N Ryman. Estimating genetic drift and effective population size fromtemporal shifts in dominant gene marker frequencies. Molecular Ecology, 8(7):1171–1178,1999.[88] Masatoshi Nei and Fumio Tajima. Genetic drift and estimation of effective population size.Genetics, 98(3):625–640, 1981.[89] Jorge Nocedal and Stephen Wright. Numerical optimization. Springer Science & BusinessMedia, 2006.103[90] Angela Burleigh, Steven McKinney, Jazmine Brimhall, Damian Yap, Peter Eirew, StevenPoon, Viola Ng, Adrian Wan, Leah Prentice, Lois Annab, J Carl Barrett, Carlos Caldas,Connie Eaves, and Samuel Aparicio. A co-culture genome-wide RNAi screen with mammaryepithelial cells reveals transmembrane signals required for growth and differentiation. BreastCancer Res., 17:4, January 2015.[91] Alejandra Bruna, Oscar M Rueda, Wendy Greenwood, Ankita Sati Batra, Maurizio Callari,Rajbir Nath Batra, Katherine Pogrebniak, Jose Sandoval, John W Cassidy, Ana Tufegdzic-Vidakovic, Stephen-John Sammut, Linda Jones, Elena Provenzano, Richard Baird, PeterEirew, James Hadfield, Matthew Eldridge, Anne McLaren-Douglas, Andrew Barthorpe,Howard Lightfoot, Mark J O’Connor, Joe Gray, Javier Cortes, Jose Baselga, ElisabettaMarangoni, Alana L Welm, Samuel Aparicio, Violeta Serra, Mathew J Garnett, and Car-los Caldas. A biobank of breast cancer explants with preserved intra-tumor heterogeneity toscreen anticancer compounds. Cell, 167(1):260–274.e22, September 2016.[92] Peter M Haverty, Eva Lin, Jenille Tan, Yihong Yu, Billy Lam, Steve Lianoglou, Richard MNeve, Scott Martin, Jeff Settleman, Robert L Yauch, et al. Reproducible pharmacogenomicprofiling of cancer cell line panels. Nature, 533(7603):333–337, 2016.[93] Wesley Tansey, Kathy Li, Haoran Zhang, Scott W Linderman, Raul Rabadan, David M Blei,and Chris H Wiggins. Dose-response modeling in high-throughput cancer drug screenings: Acase study with recommendations for practitioners. arXiv preprint arXiv:1812.05691, 2018.[94] Montserrat Rojo de la Vega, Montserrat Rojo de la Vega, Eli Chapman, and Donna D Zhang.NRF2 and the hallmarks of cancer, 2018.[95] Cancer Genome Atlas Network. Comprehensive molecular portraits of human breast tumours.Nature, 490(7418):61–70, October 2012.[96] Ann-Marie Patch, Elizabeth L Christie, Dariush Etemadmoghadam, Dale W Garsed, JoshyGeorge, Sian Fereday, Katia Nones, Prue Cowin, Kathryn Alsop, Peter J Bailey, Karin SKassahn, Felicity Newell, Michael C J Quinn, Stephen Kazakoff, Kelly Quek, CharlotteWilhelm-Benartzi, Ed Curry, Huei San Leong, Australian Ovarian Cancer Study Group, Anne104Hamilton, Linda Mileshkin, George Au-Yeung, Catherine Kennedy, Jillian Hung, Yoke-EngChiew, Paul Harnett, Michael Friedlander, Michael Quinn, Jan Pyman, Stephen Cordner,Patricia O’Brien, Jodie Leditschke, Greg Young, Kate Strachan, Paul Waring, Walid Azar,Chris Mitchell, Nadia Traficante, Joy Hendley, Heather Thorne, Mark Shackleton, David KMiller, Gisela Mir Arnau, Richard W Tothill, Timothy P Holloway, Timothy Semple, IvonHarliwong, Craig Nourse, Ehsan Nourbakhsh, Suzanne Manning, Senel Idrisoglu, TimothyJ C Bruxner, Angelika N Christ, Barsha Poudel, Oliver Holmes, Matthew Anderson, ConradLeonard, Andrew Lonie, Nathan Hall, Scott Wood, Darrin F Taylor, Qinying Xu, J LynnFink, Nick Waddell, Ronny Drapkin, Euan Stronach, Hani Gabra, Robert Brown, AndreaJewell, Shivashankar H Nagaraj, Emma Markham, Peter J Wilson, Jason Ellul, Orla Mc-Nally, Maria A Doyle, Ravikiran Vedururu, Collin Stewart, Ernst Lengyel, John V Pearson,Nicola Waddell, Anna deFazio, Sean M Grimmond, and David D L Bowtell. Whole-genomecharacterization of chemoresistant ovarian cancer. Nature, 521(7553):489–494, May 2015.[97] Yilong Li, Nicola D Roberts, Jeremiah A Wala, Ofer Shapira, Steven E Schumacher, KiranKumar, Ekta Khurana, Sebastian Waszak, Jan O Korbel, James E Haber, Marcin Imielinski,PCAWG Structural Variation Working Group, Joachim Weischenfeldt, Rameen Beroukhim,Peter J Campbell, and PCAWG Consortium. Patterns of somatic structural variation inhuman cancer genomes. Nature, 578(7793):112–121, February 2020.[98] Serena Nik-Zainal, Helen Davies, Johan Staaf, Manasa Ramakrishna, Dominik Glodzik, Xue-qing Zou, Inigo Martincorena, Ludmil B Alexandrov, Sancha Martin, David C Wedge, PeterVan Loo, Young Seok Ju, Marcel Smid, Arie B Brinkman, Sandro Morganella, Miriam RAure, Ole Christian Lingjærde, Anita Langerød, Markus Ringne´r, Sung-Min Ahn, SandrineBoyault, Jane E Brock, Annegien Broeks, Adam Butler, Christine Desmedt, Luc Dirix, SergeDronov, Aquila Fatima, John A Foekens, Moritz Gerstung, Gerrit K J Hooijer, Se Jin Jang,David R Jones, Hyung-Yong Kim, Tari A King, Savitri Krishnamurthy, Hee Jin Lee, Jeong-Yeon Lee, Yilong Li, Stuart McLaren, Andrew Menzies, Ville Mustonen, Sarah O’Meara, IrisPauporte´, Xavier Pivot, Colin A Purdie, Keiran Raine, Kamna Ramakrishnan, F Germa´nRodr´ıguez-Gonza´lez, Gilles Romieu, Anieta M Sieuwerts, Peter T Simpson, Rebecca Shep-herd, Lucy Stebbings, Olafur A Stefansson, Jon Teague, Stefania Tommasi, Isabelle Treilleux,105Gert G Van den Eynden, Peter Vermeulen, Anne Vincent-Salomon, Lucy Yates, Carlos Cal-das, Laura Van’t Veer, Andrew Tutt, Stian Knappskog, Benita Kiat Tee Tan, Jos Jonkers,A˚ke Borg, Naoto T Ueno, Christos Sotiriou, Alain Viari, P Andrew Futreal, Peter J Campbell,Paul N Span, Steven Van Laere, Sunil R Lakhani, Jorunn E Eyfjord, Alastair M Thomp-son, Ewan Birney, Hendrik G Stunnenberg, Marc J van de Vijver, John W M Martens,Anne-Lise Børresen-Dale, Andrea L Richardson, Gu Kong, Gilles Thomas, and Michael RStratton. Landscape of somatic mutations in 560 breast cancer whole-genome sequences.Nature, 534:47, May 2016.[99] C Curtis, S P Shah, S F Chin, G Turashvili, O M Rueda, M J Dunning, D Speed, A G Lynch,S Samarajiwa, Y Yuan, S Gra¨f, G Ha, G Haffari, A Bashashati, R Russell, S McKinney,METABRIC Group, A Langerod, A Green, E Provenzano, G Wishart, S Pinder, P Watson,F Markowetz, L Murphy, I Ellis, A Purushotham, A L Borresen-Dale, J D Brenton, S Tavare´,C Caldas, and S Aparicio. The genomic and transcriptomic architecture of 2,000 breasttumours reveals novel subgroups. Nature, 486(7403):346–352, June 2012.[100] ICGC/TCGA Pan-Cancer Analysis of Whole Genomes Consortium. Pan-cancer analysis ofwhole genomes. Nature, 578(7793):82–93, February 2020.[101] Mieke R Van Bockstal, Marie Colombe Agahozo, Ronald van Marion, Peggy N Atmodimedjo,Hein FBM Sleddens, Winand NM Dinjens, Lindy L Visser, Esther H Lips, Jelle Wesseling,and Carolien HM van Deurzen. Somatic mutations and copy number variations in breastcancers with heterogeneous HER2 amplification. Molecular oncology, 14(4):671–685, 2020.[102] Zaibing Li, Xiao Zhang, Chenxin Hou, Yuqing Zhou, Junli Chen, Haoyang Cai, Yifeng Ye,Jinping Liu, and Ning Huang. Comprehensive identification and characterization of somaticcopy number alterations in triple-negative breast cancer. International journal of oncology,56(2):522–530, 02 2020.[103] Arne V. Pladsen, Gro Nilsen, Oscar M. Rueda, Miriam R. Aure, Ørnulf Borgan, Knut Liestøl,Valeria Vitelli, Arnoldo Frigessi, Anita Langerød, Anthony Mathelier, Tone F. Bathen, ElinBorgen, Anne-Lise Børresen-Dale, Olav Engebr˚aten, Britt Fritzman, Øystein Garred, Ju¨rgenGeisler, Gry Aarum Geitvik, Solveig Hofvind, Vessela Kristensen, Rolf K˚aresen, Ole Chris-106tian Lingjærde, Gunhild Mari Mælandsmo, Bjørn Naume, Hege G. Russnes, Kristine KleiviSahlberg, Torill Sauer, Helle Kristine Skjerven, Ellen Schlichting, Therese Sørlie, David C.Wedge, Peter Van Loo, Carlos Caldas, Hege G. Russnes, and OSBREAC. DNA copy numbermotifs are strong and independent predictors of survival in breast cancer. CommunicationsBiology, 3(1):153, 2020.[104] Xin Shao, Ning Lv, Jie Liao, Jinbo Long, Rui Xue, Ni Ai, Donghang Xu, and Xiaohui Fan.Copy number variation is highly correlated with differential gene expression: a pan-cancerstudy. BMC Medical Genetics, 20(1):175, 2019.[105] Craig M Bielski, Ahmet Zehir, Alexander V Penson, Mark T A Donoghue, Walid Chatila,Joshua Armenia, Matthew T Chang, Alison M Schram, Philip Jonsson, Chaitanya Bandla-mudi, Pedram Razavi, Gopa Iyer, Mark E Robson, Zsofia K Stadler, Nikolaus Schultz, JoseBaselga, David B Solit, David M Hyman, Michael F Berger, and Barry S Taylor. Genome dou-bling shapes the evolution and prognosis of advanced cancers. Nat. Genet., 50(8):1189–1195,August 2018.[106] Samuel F Bakhoum, Bryan Ngo, Ashley M Laughney, Julie-Ann Cavallo, Charles J Mur-phy, Peter Ly, Pragya Shah, Roshan K Sriram, Thomas B K Watkins, Neil K Taunk, Mer-cedes Duran, Chantal Pauli, Christine Shaw, Kalyani Chadalavada, Vinagolu K Rajasekhar,Giulio Genovese, Subramanian Venkatesan, Nicolai J Birkbak, Nicholas McGranahan, MarkLundquist, Quincey LaPlant, John H Healey, Olivier Elemento, Christine H Chung, Nancy YLee, Marcin Imielenski, Gouri Nanjangud, Dana Pe’er, Don W Cleveland, Simon N Powell,Jan Lammerding, Charles Swanton, and Lewis C Cantley. Chromosomal instability drivesmetastasis through a cytosolic DNA response. Nature, 553(7689):467–472, January 2018.[107] Teresa Davoli, Hajime Uno, Eric C Wooten, and Stephen J Elledge. Tumor aneuploidycorrelates with markers of immune evasion and with reduced response to immunotherapy.Science, 355(6322), January 2017.[108] Neil T Umbreit, Cheng-Zhong Zhang, Luke D Lynch, Logan J Blaine, Anna M Cheng, RichardTourdot, Lili Sun, Hannah F Almubarak, Kim Judge, Thomas J Mitchell, Alexander Spektor,107and David Pellman. Mechanisms generating cancer genome complexity from a single celldivision error. Science, 368(6488), April 2020.[109] Ruli Gao, Alexander Davis, Thomas O McDonald, Emi Sei, Xiuqing Shi, Yong Wang, Pei-Ching Tsai, Anna Casasent, Jill Waters, Hong Zhang, Funda Meric-Bernstam, FranziskaMichor, and Nicholas E Navin. Punctuated copy number evolution and clonal stasis in triple-negative breast cancer. Nat. Genet., 48(10):1119–1130, October 2016.[110] Ahmet Acar, Daniel Nichol, Javier Fernandez-Mateos, George D Cresswell, Iros Barozzi,Sung Pil Hong, Nicholas Trahearn, Inmaculada Spiteri, Mark Stubbs, Rosemary Burke, AdamStewart, Giulio Caravagna, Benjamin Werner, Georgios Vlachogiannis, Carlo C Maley, LucaMagnani, Nicola Valeri, Udai Banerji, and Andrea Sottoriva. Exploiting evolutionary steeringto induce collateral drug sensitivity in cancer, 2020.[111] Kathleen Sprouffske, Lauren MF Merlo, Philip J Gerrish, Carlo C Maley, and Paul DSniegowski. Cancer in light of experimental evolution. Current Biology, 22(17):R762–R771,2012.[112] Gregory I Lang, Daniel P Rice, Mark J Hickman, Erica Sodergren, George M Weinstock,David Botstein, and Michael M Desai. Pervasive genetic hitchhiking and clonal interferencein forty evolving yeast populations. Nature, 500(7464):571–574, 2013.[113] Ketevan Chkhaidze, Timon Heide, Benjamin Werner, Marc J Williams, Weini Huang, GiulioCaravagna, Trevor A Graham, and Andrea Sottoriva. Spatially constrained tumour growthaffects the patterns of clonal selection and neutral drift in cancer genomic data. PLoS com-putational biology, 15(7):e1007243, 2019.[114] Allen W Zhang, Andrew McPherson, Katy Milne, David R Kroeger, Phineas T Hamilton,Alex Miranda, Tyler Funnell, Nicole Little, Camila PE de Souza, Sonya Laan, et al. Interfacesof malignant and immunologic clonal dynamics in ovarian cancer. Cell, 173(7):1755–1769,2018.[115] Shang-Hung Chen and Jang-Yang Chang. New insights into mechanisms of cisplatin resis-108tance: from tumor cell to microenvironment. International journal of molecular sciences,20(17):4136, 2019.[116] Hadi T. Nia, Lance L. Munn, and Rakesh K. Jain. Physical traits of cancer. Science,370(6516), 2020.[117] Luca Cavallone, Adriana Aguilar-Mahecha, Josiane Lafleur, Susie Brousse, Mohammed Al-damry, Talia Roseshter, Cathy Lan, Najmeh Alirezaie, Eric Bareke, Jacek Majewski, Cris-tiano Ferrario, Saima Hassan, Federico Discepola, Carole Seguin, Catalin Mihalcioiu, Eliz-abeth A. Marcus, Andre´ Robidoux, Jose´e-Anne Roy, Manuela Pelmus, and Mark Basik.Prognostic and predictive value of circulating tumor DNA during neoadjuvant chemotherapyfor triple negative breast cancer. Scientific Reports, 10(1):14704, 2020.[118] Rabih Said, Nicolas Guibert, Geoffrey R Oxnard, and Apostolia M Tsimberidou. Circulatingtumor DNA analysis in the era of precision oncology. Oncotarget, 11(2):188–211, 01 2020.[119] Kaustav Bera, Kurt A Schalper, David L Rimm, Vamsidhar Velcheti, and Anant Madabhushi.Artificial intelligence in digital pathology—new tools for diagnosis and precision oncology.Nature Reviews Clinical Oncology, 16(11):703–715, 2019.[120] Francisco Azuaje, Sang-Yoon Kim, Daniel Perez Hernandez, and Gunnar Dittmar. Con-necting histopathology imaging and proteomics in kidney cancer through machine learning.bioRxiv, 2019.[121] Yu Fu, Alexander W Jung, Ramon Vin˜as Torne, Santiago Gonzalez, Harald Vo¨hringer,Mercedes Jimenez-Linan, Luiza Moore, and Moritz Gerstung. Pan-cancer computationalhistopathology reveals mutations, tumor composition and prognosis. bioRxiv, 2019.[122] Mario Zanfardino, Monica Franzese, Katia Pane, Carlo Cavaliere, Serena Monti, Giusep-pina Esposito, Marco Salvatore, and Marco Aiello. Bringing radiomics into a multi-omicsframework for a comprehensive genotype–phenotype characterization of oncological diseases.Journal of Translational Medicine, 17(1):337, 2019.[123] Tadas Baltrusˇaitis, Chaitanya Ahuja, and Louis-Philippe Morency. Multimodal machine109learning: A survey and taxonomy. IEEE Transactions on Pattern Analysis and MachineIntelligence, 41(2):423–443, 2018.[124] Yixin Wang and David M Blei. The blessings of multiple causes. Journal of the AmericanStatistical Association, (just-accepted):1–71, 2019.[125] Yoshua Bengio, Tristan Deleu, Nasim Rahaman, Rosemary Ke, Se´bastien Lachapelle, OlexaBilaniuk, Anirudh Goyal, and Christopher Pal. A meta-transfer objective for learning todisentangle causal mechanisms. arXiv preprint arXiv:1901.10912, 2019.[126] Lucy L Gao, Jacob Bien, and Daniela Witten. Are clusterings of multiple data views inde-pendent? Biostatistics, 21(4):692–708, 2020.[127] Qian Du, Saul A Bert, Nicola J Armstrong, C Elizabeth Caldon, Jenny Z Song, Shalima SNair, Cathryn M Gould, Phuc-Loi Luu, Timothy Peters, Amanda Khoury, et al. Replicationtiming and epigenome remodelling are associated with the nature of chromosomal rearrange-ments in cancer. Nature communications, 10(1):1–15, 2019.110AppendixA.1 TNBC PDX experimental design diagrams.Figure A.1: The experimental design diagram for TNBC-SA609 PDX. Each node represents amouse/timepoint and edges from left to right denote parent-child relationships where tumour ma-terial from a parent node is extracted and transplanted onto the child node. The 4-digit numberon a node describes its sample ID or barcode. Nodes with a green background are treated withCisplatin while all others are left untreated. The nodes border colour denotes the type of sequenc-ing assay available. All mice designated X4 or later share ancestry with node 1584. This graphrepresents a tree. The nodes on the path from 0290 to 2454 comprise the 10 timepoints in theoriginal TNBC-SA609 branch. Note that 1584 to 1721 is established from fresh material while allother branches at timepoint X4 are derived from frozen material at least 1 calendar year later. Theoriginal timeseries was111Figure A.2: The experimental design diagram for TNBC-SA535 PDX. The design is similar toFigure A.1. Note that one major untreated branch (2895 to 3776) and one major treated branchare established 2895 to 3696.112Figure A.3: The experimental design diagram for TNBC-SA1035 PDX. The design is similar toFigure A.1. In Chapter 4, the untreated branch contains nodes from 2879 to 3631 while the treatedbranch comprises nodes 3015 to 3425.A.2 FiguresFigure A.4: the distribution of mismatch rate defined as the fraction of cells that have a mismatchbetween the inferred and jitter-fixed value of a marker.113−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−SA1035−RxSA535−Un SA535−Rx SA1035−UnSA609 mixture b SA609−b−Un SA609−RxHER2+ SA609−a−Un SA609 mixture ap53 wildtype p53 −/− a p53 −/− b0.50 0.75 1.00 1.25 1.50 1.750.50 0.75 1.00 1.25 1.50 1.75 0.50 0.75 1.00 1.25 1.50 1.75ABCDEFGHIJKABCDEFGHIJKABCDEFGHIJKABCDEFGHIJKABCDEFGHIJKCloneFigure A.5: Credible interval for the selection coefficients of all clones across all datasets. Thedark orange colour shows the posterior median value. The narrow grey line shows the 95% highestposterior density interval (HDI) while the thicker dark grey line shows the 50% HDI.114CABDEC A B D Eji0.000.250.500.751.00P(si > sj)p53 wildtypeHIEDBCFJGAH I E D B C F J G Aji0.000.250.500.751.00P(si > sj)p53 −/− aGCHFABDEG C H F A B D Eji0.000.250.500.751.00P(si > sj)p53 −/− bBCDB C Dji0.000.250.500.751.00P(si > sj)HER2+BDGFHAEB D G F H A Eji0.000.250.500.751.00P(si > sj)SA609−a−UnBDGHEFAB D G H E F Aji0.000.250.500.751.00P(si > sj)SA609 mixture aDFHGED F H G Eji0.000.250.500.751.00P(si > sj)SA609 mixture bAGDBHA G D B Hji0.000.250.500.751.00P(si > sj)SA609−b−UnDEGHFBAD E G H F B Aji0.000.250.500.751.00P(si > sj)SA609−RxAECGDFA E C G D Fji0.000.250.500.751.00P(si > sj)SA535−UnGCEFDAG C E F D Aji0.000.250.500.751.00P(si > sj)SA535−RxFGJIKHCDBEF G J I K H C D B Eji0.000.250.500.751.00P(si > sj)SA1035−UnJFEDBCKIGHJ F E D B C K I G Hji0.000.250.500.751.00P(si > sj)SA1035−RxFigure A.6: Posterior order matrix for selective coefficients of all clones in all datasets. Pi,j =P (si ≤ sj) shows the posterior probability that clone i has higher selective coefficient than clone j,with the stronger purple hues (close to 1.0) representing a higher confidence that clone i dominatesclone j, and conversely the stronger grey hues (close to 0.0) denote that clone j dominates clone i.Colours closer to white (0.5) represent no dominance. Note that for the lower diagonal elementsP (sj ≤ si) = 1− P (si ≤ sj) and are omitted for clarity. The diagonal entries are to guide the eyesonly.115Figure A.7: Overview of experimental design and PDX growth curves. A) Top: Serially passaged184-hTERT L9 WT cell line; Bottom: 184-hTERT L9 95.22; p53 null cell line by CRISPR tech-nology. Parallel branches P53−/− a and P53−/− b were derived from the same tenth passage. B)Western blot confirming knock out of TP53 from 184-hTERT WT cell line. Clones shown alongtop. C) Top: Schematic for PDX timeseries; Bottom: Serial sampling of HER2+ and TNBC PDXtumours; Dark grey circles represent each sampled mouse for scWGS. The light grey circles repre-sent the replicates of tumour-bearing mice at the same timepoint. D) Individual tumour growthfrom each passage of TNBC and HER2+ PDXs. E) IHC of HER2+ and TNBC tumours at earlyand late passages, 4x and 20x (insets). Scale bars 500 µm and 100 µm (insets). Antibodies andTMA scores.116Figure A.8: Single nucleotide variant (SNV) and Copy number profiles of the cell lines and PDXlibraries. A) Phylogenetic tree with clone labels next to SNV and copy number profile heatmapsfor p53-/-a. Rows are individual cells and columns are SNVs (left) and copy number alterations(right). For SNVs, orange indicates high probability while blue indicates low probability. For copynumber, legend encoding copy number colours is shown (right). B-E) Similar to A) but for p53-/-b,HER2+, TNBC, and TNBC Rx respectively117Figure A.9: TNBC-SA609 PDX clonal dynamics with and without treatment. a) Heatmap representation of copy number profiles of 841cells, grouped in 6 phylogenetic clades. b) Phylogeny as per Figure 4.2 for TNBC-SA609 PDX untreated branch. c) Observed clonalabundances, and d) distribution over magnitude of difference between selective coefficients of pairs of clones. e-h) Analogous plots forthe treated branch (n=1,593 cells).118Figure A.10: Structural variant and mutation rates of PDX lines. Distribution over copy numberbreakpoints/cell as a function of generation for A) TNBC B) HER2+. Clone specific distributionsover copy number breakpoints/cell, coloured by fitness coefficients for C) TNBC D) HER2+. Clonespecific distributions over point mutations/cell, coloured by fitness coefficients for E) HER2+ andF) TNBC119Figure A.11: Fitness validation with tumour mixing and drug perturbation. A) Schematic overviewof clonal mixture experiment showing source samples from the original timeseries and serial prop-agation into a new line. B) Mouse body weight graph recorded during maximum tolerated dose(MTD) evaluation of cisplatin in NRG mice (n=3 in each study cohort). C) Experimental design ofcisplatin treatment in PDX.The residual tumour from one treated mouse was re-transplanted in thenext (n=4). The solid blue colour represent cisplatin treated tumours (UT,UTT,UTTT,UTTTT);blue outlined in grey represents drug holiday (UTU,UTTU,UTTTU). Grey represent the untreatedseries (U,UU,UUU,UUUU,UUUUU). D) Tumour response curves in each cycle of cisplatin treat-ment. 120A.3 TablesDataset Timepoint1 Timepoint2 Normalised counts Copy number rate1 SA501 X2 X5 0.02 0.012 SA501 X5 X6 0.02 0.023 SA501 X6 X8 0.19 0.094 SA501 X8 X11 0.05 0.025 SA501 X11 X15 0.04 0.016 SA535 X1 X5 0.17 0.047 SA535 X5 X8 0.34 0.11Table A.2: Estimated CNA rate over time in timeseries datasets used in Chapter 2. Note that theCopy number rate is not normalised by the length of the passage (in number of generations) ormean ploidy.Dataset Ne Number of SNVs1 p53 WT 985.79 NA2 p53-/-a 614.93 6,6453 p53-/-b 422.14 7,9054 HER2+ 1,461.70 40,3095 TNBC 468.30 33,1786 TNBC-Mixture 177.01 NA7 TNBC-Rx 333.34 25,685Table A.3: Parameters per dataset. Effective population size estimates and number of loci used inthe SNV analysis.Dataset  runtime niter ∆τ avg ESS min ESS rejection rate1 p53 WT 0.01 2:57:22. 10,000 0.05 911.90 448.19 0.772 p53-/-a 0.01 35:18:46. 100,000 0.05 3209.54 506.66 0.613 p53-/-b 0.03 31:29:55. 100,000 0.04 2077.27 906.95 0.894 HER2+ 0.01 48:57:27. 100,000 0.05 26303.53 25299.69 0.385 TNBC 0.03 40:05:32. 100,000 0.05 1379.81 223.08 0.916 TNBC-mixture 0.01 29:13:01. 100,000 0.04 3745.76 1268.09 0.737 TNBC-Rx 0.03 30:05:33. 100,000 0.10 4362.86 2768.55 0.738 WGS-bulk-TNBC 0.012 35:12:29. 100,000 0.10 14560.09 11754.76 0.71Table A.4: Real-world data parameters.121dataset sample.id lib.id tp total !mouse qual. !sphase !lmr final finalr1 SA501 SA501X2XB00096 A95621B X2 623 570 524 580 568 498 4982 SA501 SA501X2XB00097 A96171A X2 492 36 21 474 37 20 203 SA501 SA501X2XB00097 A96109A X2 656 59 37 399 96 34 344 SA501 SA501X5XB00877 A95670A X5 615 326 253 575 252 238 2285 SA501 SA501X6XB00294 A96213A X6 1000 326 152 886 168 62 626 SA501 SA501X6XB00969 A95670B X6 636 383 325 602 300 319 2727 SA501 SA501X8XB01694 A90696ABC X8 3672 3517 273 3482 1985 261 2418 SA501 SA501X11XB00529 A96187A X11 1063 1024 878 888 958 779 7639 SA501 SA501X11XB00529 A96174A X11 1234 1192 291 1007 258 116 5810 SA501 SA501X15XB00929 A96173A X15 1388 997 558 1131 667 399 39911 TNBC-SA535 SA535X1XB00174 A96165B X1 379 208 73 218 277 62 6212 TNBC-SA535 SA535X5XB00517 A95732A X5 928 613 372 776 464 332 32413 TNBC-SA535 SA535X8XB00143 A95736A X8 1072 410 194 956 354 154 15214 OVA-2295 SA921 A90554A TOV2295(R) 628 506 415 524 525 392 39215 OVA-2295 SA922 A90554B OV2295(R2) 596 466 303 518 467 296 29616 OVA-2295 SA1090 A96213A OV2295 1484 1438 1144 1092 1402 850 838Table A.5: Summary of real-world datasets used. final is the final number of cells after all filters except for !lmr are applied. finaladditionally filters out lmr cells, those that have total mapped reads fewer than 500,000. Abbreviations used are tp: time point; qual. :quality; !sphase: not S-phase; lmr: low mapped reads; !lmr: not low mapped reads.122A.4 Preprocessing of the CNA data for phylogenetic inferenceGet copy number valuesCNA states are stored in cn.csv and this file is the input to our preprocessing pipeline. For thiswork the data was stored in the cloud (Microsoft Azure) and the scgenome application programminginterface (API) suite was used to access and download the data. Please see https://github.com/shahcompbio/scgenome for documentation. This API was developed by authors in reference [20].The scgenome API ensures that only cells with the correct sample ids are selected, and removescontrol cells and cells that have fewer than 10,000 mapped reads. The API also provides for eachcell a quality score that we call SCG-score. This score ranges from zero to one where values closerto one indicate higher quality cells. Lower quality cells are expected to be noisy, that is to havemany non-integer CNA values and very low read counts. See the [20] for more details.Drop low-mappability binsSome copy number bins are located at parts of the genome where sequencing is difficult, for exampledue to inaccessibility of the genome at that position. These bins are sometimes assigned imputedCN states that represent a sequencing artefact. This is reflected in their mappability score. We filterthe CNA matrix to keep high-mapability bins cn bin filtered.csv. In this work we use a cutoffthreshold of map >=.99 that yields 4375/6206 or 70.5% of the bins. The list of kept bins is identicalacross all datasets. Low-mapbaility bins may constitute phylogenetic markers that are not reflectiveof the underlying biology and result in grouping of cells that share these sequencing artefacts but arenot evolutionarily related. Moreover, using the 0.99 cutoff, we still recover hundreds of phylogeneticmarkers that are present in at least 5% of cells. A potential pitfall can occur in datasets that havevery few copy number changepoints to accurately distinguish subpopulations. In such cases filteringtoo many genomic bins may result in missing some subpopulations. A lack of ground-truth datasetsmakes it difficult to set this threshold in a principled way.Drop low-quality cellsIn this step, a second round of quality control is done. Cells with SCG-scores of over 0.75 arekept [20], those that are suspected to be contaminated (e.g., mouse cells) or cycling cells are123removed. This results in cn bin cell filtered.csv. The 0.75 cutoff is chosen as it balancessensitivity and specificity of the single cell scoring algorithm.Drop cells with excess CNA changesSome cells show a jumpy CNA profile in which there are excessive copy number changes. It appearsthat these cells are either in early or late stages of division and were missed by the scgenome API.The CNA profile of early replicating cells is patterned by seemingly scattered focal amplificationswhile the late replicating cells show scattered focal deletions. Note that not all parts of the genomeduplicate at the same time during mitosis [127]. Regions that start duplicating later will show ashaving focal deletions in cells captured at their later replicating stage; these regions would not havestarted to duplicate by the time the sample was prepared for sequencing. One possible explanationis how the read counts in cells are normalised to integer CN states in the DLP+ platform. Whenestimating the CN states of cells, the normalizing procedure assumes that the most likely ploidyfor all cells is diploid [20]. For each cell, the observed read counts are rescaled as closely to thediploid state as possible, while retaining integer CN states. These scattered patterns (Figure A.12)are hypothesised to not directly reflect the evolutionary history of the cells and are detrimental tophylogenetic tree inference.2.5. Tree evaluationFigure 2.8: Phylogenetic tree and CNA profile heatmap for the SA501 dataset.2.5 Tree evaluation2.5.1 Predictive testTo evaluate the inferred trees, we suggest a test that involves predicting the entries in the inputbinary matrix given to the tree inference method. We take the binarised input matrix yC⇥L, theinput matrix to the sitka algorithm as described in Section 2.3.4 as ground truth. Consider aninferred tree, ⌧ , and the corresponding genotype matrix g = T (⌧). In general the inferred treesfrom the baseline methods do not have named internal nodes, nor do they have the same number ofinternal nodes as the number of loci L. Therefore we do not know which locus in the inferred tree⌧ corresponds to which locus in the matrix y. We note that this is not the case with trees inferredfrom sitka where the internal nodes of the tree correspond to the columns of the induced genotypematrix. As a result, for methods other than sitka, for each column in the input data matrix, wepick a clade in ⌧ that has the highest prediction accuracy for the entries in that column.For each method, we report Youden’s J index [109] which is equal to sensitivity + specificity -1. We define below the function h to be a binary classification counts matrix, i.e., for two C-vectorsw and z it forms the confusion matrix. h : {0, 1}1⇥C ⇥ {0, 1}1⇥C 7! {0, 1}2⇥2 wherehi,j(w, z) =Xc2C1 (wc = i)1 (zc = j) .25Copy number stateSporadic localised deletions on chromosome 4 (top) and chromosomes 7 and 8 (bottom) are likely a sign of late-replicating cells. CellsChromosomesClade comprising likely late-replicating cells Figure A.12: An example of replicating cells. The CNA matrix and the consensus tree show cellsfrom a HER2+ PDX dataset. Cells with excess CNA changes are not filtered out (see Section A.4).Note that a subset of cells, corresponding to the green clade on the tree, exhibit scattered localiseddeletions (the arrows on the right) across the genome. These late replicating cells form a finger -likeclade in the tree. The similarity in the CNA profiles is not due to common evolutionary historyand th r fore they should be removed (see Section A.4). The top inset shows chromosome 4 whilethe bottom inset spans chromosomes 7 and 8.124Here we rank cells by the number of changes in their copy number states (a change is measuredbetween consecutive bins) and pick the bottom 90-th percentile. This threshold was set via visualinspection. The file cn bin cell filtered no jump.csv contains the integer copy number statewith the final list of cells and genomic bins. An example input matrix is shown in Figure 2.4-awhere the integer copy numbers are colour-coded in a heatmap. Attrition rate due to filtering ofcells is shown in Table A.5.A.5 Computing CNA rate over timeComputing the CNA rate is non-trivial. The reference timepoint is heterogenous and therefore toavoid conflating selection with mutation rate, we pick a set of markers that are not present in thereference timepoint (cellular prevalence <0.001) and call it L0. In each timepoint i (other than thereference timepoint), find the subset of markers in L0 that are present in that timepoint (e.g., havea cellular prevalence >0.01) and call this subset Li. Now for each pair of consecutive timepoints,count the number of markers that are not identical (i.e., Li \ Li+1), normalised by the number ofmarkers in each dataset. Table A.2 shows the normalised CNA rates for the timeseries datasetsused in Chapter 2.125
Thesis/Dissertation
eng
Bioinformatics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Quantitative fitness modelling in cancer using single cell timeseries population dynamics
Text