Science, Faculty of
Computer Science, Department of
DSpace
UBCV
Dorri, Fatemeh
Doctor of Philosophy - PhD
University of British Columbia
Next generation sequencing (NGS) technology provides researchers an opportunity to study cancer genomes at different resolutions. In particular, detection and interpretation of the smallest somatic changes of the genome (single nucleotide variants) are now tractable at scale. However, significant challenges in the analysis of both bulk tumour and single cell sequencing methods remain to fully exploit the advance in technology development. Two emerging areas of applying sequencing technology to better ascertain properties of cancer evolution are (i) sequencing multiple tumour biopsies from the same patient, and (ii) single cell genome sequencing. Both of these advances represent computational challenges that I address through development of novel methods in this thesis. The first proposed method (Chapter 2) incorporates prior clonal information to improve the accuracy of detecting SNVs across the genome of multiple bulk tumour samples. The second proposed method (Chapter 3) is a statistical model that exploits the underlying phylogeny of individually sequenced cells to detect SNVs in every individual cell. The latter method identifies clone specific SNVs without the requirement of deconvolving the results from bulk sequencing data. The resultant accurate detection of SNVs (Chapter 4) helps enhance insight on the evolutionary process of tumours and genetic pathways. Together, the methods provide a toolbox for comprehensive profiling of SNVs for the study of tumour dynamics.
Somatic mutation analysis for the study of clonal evolutionin cancerbyFatemeh DorriMSc, University of Waterloo, 2012A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Computer Science)The University of British Columbia(Vancouver)April 2020c© Fatemeh Dorri, 2020The following individuals certify that they have read, and recommend to the Fac-ulty of Graduate and Postdoctoral Studies for acceptance, the thesis entitled:Somatic mutation analysis for the study of clonal evolution in cancersubmitted by Fatemeh Dorri in partial fulfillment of the requirements for the de-gree of Doctor of Philosophy in Computer Science.Examining Committee:Dr. Sohrab Shah, Department of Pathology and Laboratory Medicine, Departmentof Molecular OncologyCo-SupervisorDr. Anne Condon, Department of Computer ScienceCo-SupervisorDr. Alexandre Bouchard-Côté, Department of StatisticsSupervisory Committee MemberDr. Gabriela V. Cohen Freue, Department of StatisticsUniversity ExaminerDr. Sara Mostafavi, Department of StatisticsUniversity ExamineriiAbstractNext generation sequencing (NGS) technology provides researchers an opportu-nity to study cancer genomes at different resolutions. In particular, detection andinterpretation of the smallest somatic changes of the genome (single nucleotidevariants) are now tractable at scale. However, significant challenges in the analysisof both bulk tumour and single cell sequencing methods remain to fully exploitthe advance in technology development. Two emerging areas of applying sequenc-ing technology to better ascertain properties of cancer evolution are (i) sequenc-ing multiple tumour biopsies from the same patient, and (ii) single cell genomesequencing. Both of these advances represent computational challenges that I ad-dress through development of novel methods in this thesis. The first proposedmethod (Chapter 2) incorporates prior clonal information to improve the accuracyof detecting SNVs across the genome of multiple bulk tumour samples. The sec-ond proposed method (Chapter 3) is a statistical model that exploits the underlyingphylogeny of individually sequenced cells to detect SNVs in every individual cell.The latter method identifies clone specific SNVs without the requirement of decon-volving the results from bulk sequencing data. The resultant accurate detection ofSNVs (Chapter 4) helps enhance insight on the evolutionary process of tumoursand genetic pathways. Together, the methods provide a toolbox for comprehensiveprofiling of SNVs for the study of tumour dynamics.iiiLay SummaryThe accumulation of genetic mutations disrupts the regular activity of cells andleads to development of tumours. As cancerous cells divide, their acquired mu-tations are passed down to their descendants. This thesis proposes methods thatexploit this information for detecting mutations with a resolution down to one cell.We expect the methods to be applied for the study of tumours clonal dynamics. De-ciphering the clonal dynamics of tumours can potentially lead to enhanced cancerdiagnosis and treatment.ivPrefaceChapter 2 is a modified version of the material published in Dorri et al. [Dorri,F., Jewell, S., Bouchard-Côté, A., Sohrab P. Shah. Somatic mutation detection andclassification through probabilistic integration of clonal population information.Nature Communications Biology 2, 44 (2019)]. The MuClone model was jointlyconceived by myself, Dr. Shah, and Dr. Bouchard-Côté. I developed the MuClonesoftware with Sean Jewell. I performed all computational analyses and I am theoriginal creator and copyright holder of all the figures presented in this chapter. Ico-wrote the text with Sean Jewell, Dr. Shah and Dr. Bouchard-Côté.Chapter 3 contains unpublished material that is being prepared for a submissionto a peer review journal. The model was jointly conceived by myself, Dr. Shahand Dr. Bouchard-Côté. I developed the CellMutScope software and performedall computational analyses described. The software was implemented within theCorrupt package. I wrote the text and I am the original creator and copyrightholder of all the figures presented in this chapter.Chapter 4 contains unpublished material that is being prepared for submissionto a peer reviewed journal. I performed all computational analyses described inchapter 4. I wrote the text. The figures are generated from a modified version ofTyler Funnel’s code for similar figures. I am the original creator and copyrightholder of all figures presented in this chapter. The data presented in this chapterfrom a manuscript in preparation with Sohrab Salehi. The xenografting data wascollected with informed patient consent according to procedures approved by theEthics Committee at the University of British Columbia, under protocols H06-00289 BCCA-TTR-BREAST and H11-01887 Neoadjuvant Xenograft Study.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xGlossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Clonal evolution in cancer and intra-tumour heterogeneity . . . . 21.2 The landscape of genetic variants in tumours . . . . . . . . . . . . 51.2.1 Single nucleotide variants . . . . . . . . . . . . . . . . . 51.2.2 Structural variants . . . . . . . . . . . . . . . . . . . . . 61.2.3 SNVs and SVs in the study of tumour dynamics . . . . . . 61.3 Next generation sequencing technologies . . . . . . . . . . . . . 61.3.1 Bulk genome sequencing . . . . . . . . . . . . . . . . . . 9vi1.3.2 Single cell sequencing . . . . . . . . . . . . . . . . . . . 101.4 Literature review on SNV calling models . . . . . . . . . . . . . . 111.4.1 Literature review on SNV calling from bulk sequencing data 111.4.2 Literature review on SNV detection from single cell se-quencing data . . . . . . . . . . . . . . . . . . . . . . . . 151.5 Motivations and research contribution . . . . . . . . . . . . . . . 162 Somatic mutation detection and classification through probabilisticintegration of clonal population structure . . . . . . . . . . . . . . . 182.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.2 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . 202.3 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.3.1 Problem formulation . . . . . . . . . . . . . . . . . . . . 212.3.2 Model description . . . . . . . . . . . . . . . . . . . . . 212.3.3 Inference . . . . . . . . . . . . . . . . . . . . . . . . . . 252.4 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . 262.4.1 Synthetic data . . . . . . . . . . . . . . . . . . . . . . . . 262.4.2 Real data . . . . . . . . . . . . . . . . . . . . . . . . . . 332.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433 Single cell somatic mutation detection through incorporation of theunderlying phylogeny . . . . . . . . . . . . . . . . . . . . . . . . . . 553.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 563.2 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . 583.3 Corrupt model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 603.3.1 Probability model . . . . . . . . . . . . . . . . . . . . . . 623.3.2 Inference . . . . . . . . . . . . . . . . . . . . . . . . . . 623.3.3 Approximation of the posterior distribution . . . . . . . . 643.3.4 Summary via minimum Bayes estimator . . . . . . . . . . 653.4 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 653.4.1 Incorporating single nucleotide data in the Corrupt model 653.4.2 Detection of SNVs at every individual cell . . . . . . . . . 683.4.3 Model’s asymptotic complexity . . . . . . . . . . . . . . 72vii3.5 Benchmarking experiments . . . . . . . . . . . . . . . . . . . . . 723.5.1 Synthetic data generation . . . . . . . . . . . . . . . . . . 723.5.2 Synthetic data evaluation . . . . . . . . . . . . . . . . . . 733.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 774 Real data application . . . . . . . . . . . . . . . . . . . . . . . . . . 814.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 814.2 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 824.2.1 Ovarian cancer . . . . . . . . . . . . . . . . . . . . . . . 824.2.2 Breast cancer xenograft samples . . . . . . . . . . . . . . 834.3 Data analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 844.3.1 Data cleaning . . . . . . . . . . . . . . . . . . . . . . . . 844.3.2 Main analysis . . . . . . . . . . . . . . . . . . . . . . . . 844.4 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . 914.5 Mutations in high impact genes . . . . . . . . . . . . . . . . . . . 1224.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1235 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1255.1 Summary of contributions . . . . . . . . . . . . . . . . . . . . . 1255.2 Future work and discussion . . . . . . . . . . . . . . . . . . . . . 126Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129viiiList of TablesTable 2.1 Summary of high grade serous ovarian cancer data set . . . . . 36Table 2.2 Summary of the non-small cell lung cancer data set . . . . . . 40Table 2.3 Total number of false negative calls across multiple samples ofnon-small cell lung cancer patients for different algorithms. . . 41Table 4.1 Raw data information for OV 2295, SA532 and SA609 . . . . . 94Table 4.2 The total number of cells before and after quality control step . 106Table 4.3 High impact and non-synonymous coding genes in OV 2295 . . 123Table 4.4 High impact and non-synonymous coding genes in SA532 . . . 123Table 4.5 High impact and non-synonymous coding genes in SA609 . . . 124ixList of FiguresFigure 1.1 A schematic view of bulk and single cell genome sequencing . 5Figure 1.2 Workflow of next generation sequencing technology . . . . . 8Figure 2.1 The graphical model of MuClone . . . . . . . . . . . . . . . 22Figure 2.2 Distribution of the cellular prevalences across different sam-ples in multiple runs . . . . . . . . . . . . . . . . . . . . . . 27Figure 2.3 Distribution of the cellular prevalences across different sam-ples in one random run . . . . . . . . . . . . . . . . . . . . . 28Figure 2.4 A schematic view of clonal information . . . . . . . . . . . . 29Figure 2.5 MuClone’s performance with inaccurate clonal information . 30Figure 2.6 MuClone’s performance as a function of different wildtype prior 32Figure 2.7 MuClone’s performance as a function of tumour content . . . 33Figure 2.8 MuClone’s performance as a function of different depth . . . 34Figure 2.9 MuClone’s classification performance in synthetic data . . . . 35Figure 2.10 Performance comparison of different methods on whole genomesequencing data from patients with high grade serous ovariancancer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37Figure 2.11 Performance comparison of different methods on whole genomesequencing data for patient 1 with high grade serous ovariancancer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38Figure 2.12 MuClone’s Roc curves and the area under the curve (AUC) forpatient 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39xFigure 2.13 Performance comparison of different methods on whole genomesequencing data for patient 2 with high grade serous ovariancancer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40Figure 2.14 MuClone’s Roc curves and the area under the curve (AUC) forpatient 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41Figure 2.15 Performance comparison of different methods on whole genomesequencing data for patient 3 with high grade serous ovariancancer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42Figure 2.16 MuClone’s Roc curves and the area under the curve (AUC) forpatient 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43Figure 2.17 Performance comparison of different methods on whole genomesequencing data for patient 4 with high grade serous ovariancancer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44Figure 2.18 MuClone’s Roc curves and the area under the curve (AUC) forpatient 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45Figure 2.19 Performance comparison of different methods on whole genomesequencing data for patient 7 with high grade serous ovariancancer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46Figure 2.20 MuClone’s Roc curves and the area under the curve (AUC) forpatient 7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47Figure 2.21 Performance comparison of different methods on whole genomesequencing data for patient 9 with high grade serous ovariancancer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48Figure 2.22 MuClone’s Roc curves and the area under the curve (AUC) forpatient 9 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49Figure 2.23 Performance comparison of different methods on whole genomesequencing data for patient 10 with high grade serous ovariancancer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50Figure 2.24 MuClone’s Roc curves and the area under the curve (AUC) forpatient 10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51Figure 2.25 Classification of mutations of patient 1 with high grade serousovarian cancer across 6 samples . . . . . . . . . . . . . . . . 52xiFigure 2.26 Comparison of detected mutations from MultiSNV, MuTect,MuClone, and TRACERx on whole exome sequencing datafrom non-small cell lung cancer . . . . . . . . . . . . . . . . 53Figure 2.27 MuClone’s performance with inaccurate clonal information . 54Figure 3.1 A schematic representation of the phylogeny . . . . . . . . . 59Figure 3.2 A schematic view of removing and adding a node in the Cor-rupt model . . . . . . . . . . . . . . . . . . . . . . . . . . . 63Figure 3.3 A schematic representation of CNV and SNV tree . . . . . . . 71Figure 3.4 The performance (Sensitivity and Specificity) of CellMutScopeacross different data coverages . . . . . . . . . . . . . . . . . 74Figure 3.5 The performance (Youden’s index) of CellMutScope acrossdifferent data coverages . . . . . . . . . . . . . . . . . . . . 75Figure 3.6 The performance (D) of CellMutScope across different datacoverages . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75Figure 3.7 The performance (Sensitivity and Specificity) of CellMutScopeacross different number of loci used at inference step . . . . . 76Figure 3.8 The performance (Youden’s index) of CellMutScope acrossdifferent number of loci used at inference step . . . . . . . . . 76Figure 3.9 The performance (D) of CellMutScope across different num-ber of loci used at inference step . . . . . . . . . . . . . . . . 77Figure 3.10 The performance (Sensitivity and Specificity) of CellMutScopeacross different number of cells . . . . . . . . . . . . . . . . 78Figure 3.11 The performance (Youden’s index) of CellMutScope acrossdifferent number of cells . . . . . . . . . . . . . . . . . . . . 79Figure 3.12 The performance (D) of CellMutScope across different num-ber of cells . . . . . . . . . . . . . . . . . . . . . . . . . . . 80Figure 4.1 A schematic view of single cell sequencing data analysis . . . 86Figure 4.2 The distribution of the number of cells from different samplesin each clone for OV 2295. . . . . . . . . . . . . . . . . . . . 86Figure 4.3 The distribution of the number of cells from different samplesin each clone for SA532. . . . . . . . . . . . . . . . . . . . . 87xiiFigure 4.4 The distribution of the number of cells from different samplesin each clone for SA609. . . . . . . . . . . . . . . . . . . . . 87Figure 4.5 The copy number heatmap for OV 2295 data . . . . . . . . . . 88Figure 4.6 The copy number heatmap for SA532 data . . . . . . . . . . . 89Figure 4.7 The copy number heatmap for SA609 data . . . . . . . . . . . 90Figure 4.8 The SNV call probability distribution for OV 2295 . . . . . . . 92Figure 4.9 The SNV call probability distribution for SA532. . . . . . . . . 93Figure 4.10 The SNV call probability distribution for SA609 . . . . . . . . 93Figure 4.11 The SNV call probability heatmap for OV 2295 before qualitycontrol step . . . . . . . . . . . . . . . . . . . . . . . . . . . 95Figure 4.12 The number of variant reads heatmap for OV 2295 before qual-ity control step . . . . . . . . . . . . . . . . . . . . . . . . . 96Figure 4.13 The SNV call probability heatmap for SA532 before qualitycontrol step . . . . . . . . . . . . . . . . . . . . . . . . . . . 98Figure 4.14 The number of variant reads heatmap for SA532 before qualitycontrol step . . . . . . . . . . . . . . . . . . . . . . . . . . . 99Figure 4.15 The SNV call probability heatmap for SA609 before qualitycontrol step . . . . . . . . . . . . . . . . . . . . . . . . . . . 100Figure 4.16 The number of variant reads heatmap for SA609 before qualitycontrol step . . . . . . . . . . . . . . . . . . . . . . . . . . . 101Figure 4.17 The SNV call probability heatmap for OV 2295 after qualitycontrol step with threshold equals 0.001 . . . . . . . . . . . . 102Figure 4.18 The number of variant reads heatmap for OV 2295 after qualitycontrol step with threshold equals 0.001 . . . . . . . . . . . . 103Figure 4.19 The SNV call probability heatmap for OV 2295 after qualitycontrol step with threshold equals 0.005 . . . . . . . . . . . . 104Figure 4.20 The number of variant reads heatmap for OV 2295 after qualitycontrol step with threshold equals 0.005 . . . . . . . . . . . . 105Figure 4.21 The SNV call probability heatmap for SA532 after quality con-trol step with threshold equals 0.001 . . . . . . . . . . . . . . 107Figure 4.22 The number of variant reads heatmap for SA532 after qualitycontrol step with threshold equals 0.001 . . . . . . . . . . . . 108xiiiFigure 4.23 The SNV call probability heatmap for SA532 after quality con-trol step with threshold equals 0.005 . . . . . . . . . . . . . . 109Figure 4.24 The number of variant reads heatmap for SA532 after qualitycontrol step with threshold equals 0.005 . . . . . . . . . . . . 110Figure 4.25 The SNV call probability heatmap for SA609 after quality con-trol step with threshold equals 0.001 . . . . . . . . . . . . . . 111Figure 4.26 The number of variant reads heatmap for SA609 after qualitycontrol step with threshold equals 0.001 . . . . . . . . . . . . 112Figure 4.27 The SNV call probability heatmap for SA609 after quality con-trol step with threshold equals 0.005 . . . . . . . . . . . . . . 113Figure 4.28 The number of variant reads heatmap for SA609 after qualitycontrol step with threshold equals 0.005 . . . . . . . . . . . . 114Figure 4.29 The mean SNV call probability heatmap of permutation analy-sis for OV 2295 after quality control step with threshold 0.001 116Figure 4.30 The number of variant reads heatmap of permutation analysisfor OV 2295 after quality control step with threshold 0.001 . . 117Figure 4.31 The mean SNV call probability heatmap of permutation analy-sis for SA532 after quality control step with threshold 0.001 . 118Figure 4.32 The number of variant reads heatmap of permutation analysisfor SA532 after quality control step with threshold 0.001 . . . 119Figure 4.33 The mean SNV call probability heatmap of permutation analy-sis for SA609 after quality control step with threshold 0.001 . 120Figure 4.34 The number of variant reads heatmap of permutation analysisfor SA609 after quality control step with threshold 0.001 . . . 121xivGlossaryCNV Copy number variantsDLP Direct library preparationDNA Deoxyribonucleic acidFFPE Formalin-Fixed Paraffin-EmbeddedNGS Next generation sequencingPCR Polymerase chain reactionPT Parallel temperingROC Receiver operating characteristicSCS Single cell sequencingSNP Single nucleotide polymorphismSNV Single nucleotide variantsSV Structural variantsTNBC Triple-negative breast cancerVAF Variant allelic fractionWES Whole exome sequencingWGA Whole genome amplificationWGS Whole genome sequencingxvAcknowledgementsFirst and foremost, I wish to thank multitude of people who helped me throughoutmy journey as a PhD student. I would like to express my sincere gratitude to myadvisor, Prof. Sohrab Shah, for giving me an invaluable opportunity to work onchallenging and interesting projects over the past 6 years. Also, I would like tothank Prof. Anne Condon for all her support during these years. I would liketo extend my appreciation for Prof. Alexandre Bouchard-Côté for serving on mydissertation committee. It has been a pleasure to work with and learn from himduring these years.I believe that graduate school was one of the most interesting pages of my life.I would like to thank all my friends and my colleagues at BC Cancer researchCenter, Dr. Razieh Annabestani, Dr. Sayeh Rajabi, Dr. Nilgoon Zarei, and manyothers for filling this page of my life with the most wonderful memories. Thankyou all!I would like to thank all who helped with their spiritual support. I feel a deepsense of gratitude for my mother and father, Zahra and Behrouz, who formed partof my vision, and taught me things that really matter in life. Their infallible loveand support has always been my strength. Their patience and sacrifice will remainmy inspiration throughout my life. I like to thank my twin sister, Faezeh and mybrother-in-law, Hamid Mahini for their sincere help and support. They are alwaysthanked in my heart for their enormous love. I would like to thank my lovelydaughter, Elina, for all the joy and energy she gave me after her birth. Last butnot least, I am deeply grateful to the love of my life and my husband, MohsenKeshavarz Akhlaghi, for all support and sacrifice he made during every single stepof this Journey.xviIt is impossible to remember all, and I apologize to those I have inadvertentlyleft out. Lastly, thank you all and thank God!xviiTo my lovely ElinaChapter 1Introduction“It is paradoxical, yet true, to say, that the more we know, themore ignorant we become in the absolute sense, for it is only throughenlightenment that we become conscious of our limitations. Preciselyone of the most gratifying results of intellectual evolution is the con-tinuous opening up of new and greater prospects.”- Nikola TeslaCancer is an evolutionary process driven by the accumulation of genetic vari-ants. Intra-tumour heterogeneity is an inevitable consequence of this evolutionaryprocess (Section 1.1). Identification of the underlying heterogeneous populations(clones) is an important step toward understanding the constituent parts of a tu-mour. Genetic variants can occur at single nucleotides, or there may be largervariants in the genome (Section 1.2). Because single nucleotide variants (SNVs)can be used as a marker in deciphering clonal dynamics, and because they areimportant in genetic pathways, accurate SNV detection is critical. Advanced se-quencing technologies provide raw data with SNV information covered (1.3). Inspite of numerous methods developed in the field (1.4), it is still desired to use thisdata to detect SNVs with higher accuracy. Incorporating the data obtained fromdifferent sequencing technologies, we proposed methods to detect SNVs with highaccuracy, particularly SNVs that are only present in a relatively small population ofcells (Section (1.5))11.1 Clonal evolution in cancer and intra-tumourheterogeneityDuring the cell division process, the DNA content of a cell has to be accuratelytransmitted in order to maintain the genetic integrity. However, the limited fi-delity of this process results in genetic variants. Genetic variants arise in DNA ofcells through different mechanisms including DNA division process or exposureto physical or chemical agents. DNA damage is induced by both endogenous (e.g.metabolic by-products such as reactive oxygen species) and exogenous factors (e.g.UV radiation and viruses) [38]. It is instructive to note the normal haploid humangenome consists of about 3×109 base pairs. So, with the approximate 10−9 basepair mutation (single nucleotide variant) rate, it is likely that daughter cells inherita new mutation during every division process [73, 130]. Unrepaired mutations arepassed to the daughter cells and all subsequent descendants. Cells also possess aDNA repair mechanism. The rate of mutation accumulation is therefore a combinedfunction of improper DNA replication inducing mutations and defective DNA repairmechanism [15]. Defects in DNA replication process or repair mechanism, or expo-sure to exogenous factors can increase the overall rate of mutations. For example,the mutation rates respectively reach as high as 12.6, 135.1 and 173.9 mutationsper megabase of genome in high-grade serous ovarian cancer, breast cancer, andnon-small cell lung cancer [20].Most mutations (neutral or passenger mutations) do not impact a cell’s pheno-type. A vital phynotype, cell birth to death ratio, is affected by either (i) acquiringdisadvantageous mutations that slow down the division process (this causes highercell death rate), or (ii) driver mutations that cause an increase in cell birth rate.Acquiring a driver mutation can lead to a higher proliferation rate and develop-ment of cancerous cells [14, 52]. Besides genetic variants, epigenetic mechanismslike methylation and heterogeneous tumour micro-environments may also affectthe proliferation and survival of tumour cells [142].A group of cells with similar genetic variants and similar phenotypic propertiesare called clones [120]. Clonal expansion is the uncontrolled growth of a group oftumour cells with higher survival rate [48, 52]. A growing tumour cell populationincludes different (heterogeneous) clones with distinct genetic profiles [88, 92].2Intra-tumour heterogeneity describes the existence of heterogeneous clones withdistinct phenotypic properties. For example, a clone with selective advantage sur-vives better under some treatment [49, 55], or it can spread to other sites in thebody by metastasis.Intra-tumour heterogeneity is observed in various cancers including lung [27,143], breast [36, 89, 90, 108, 137], ovarian [8, 106, 123], pancreatic [16, 136], kid-ney [42, 43], colorectal [118, 122], brain [117, 121], acute lymphoblastic leukemia [4],and prostate [18, 26, 51]. Intra-tumour heterogeneity complicates the accurate di-agnosis of tumour, the prediction of disease progresses, and treatments [68, 77].For example, in non–small cell lung cancer, T790M mutation is found to be re-sistant to EGFR inhibitors [29]. KRAS mutations in colon cancers is a markerof resistance to EGFR-I therapy [112]. BRCA1 reversion mutations in breast andovarian cancers are identified as conferring resistance to PARP inhibitors and/orcisplatin [91]. ESR1 mutations are found to cause resistance in ER+ breast cancerstreated with tamoxifen [61].Identification of clones and measuring their dynamics (variations in clone’sfitness and growth) is critical in deciphering the evolutionary process of cancer.The process of cells acquiring new mutations and passing them on to the daughtercells leaves a record of ancestral relationship between cells (see Figure 1.1 panela). Basically, each tumour cell is a copy of its parents with a few probable muta-tions. Therefore, tumour evolution can be studied by inferring the genetic contentof each cell. In order to determine the genetic profile of the existing clones, tissuesamples dissected from parts of the body, are sequenced and aligned to the refer-ence genome. This provides the raw data for inferring the genetic profile and forthe subsequent studies.Single nucleotide variants are the most common type of genetic variants [66].In spite of the progress in developing SNV detection algorithms in recent years, theproblem remains challenging particularly for detecting SNVs that exist only in arelatively small subset of cells (rare clones). This thesis is focused on detectionof SNVs for the study of clonal populations and their dynamics. We showed thatinjecting clonal information improves the performance of detecting SNVs.3(a) (b)(c)(d)Figure 1.1: (Caption next page.)4Figure 1.1: (a) A hypothetical phylogenetic tree shows the ancestral relation-ship between clones (circles). Each clone is identified by mutations(stars). Unlike traditional phylogenetic trees, the internal clones maycontribute to the observed data. The diameter of the circles (clones)is proportional to the clone’s abundance. The black star represents anancestral mutation, and other stars are clone specific mutations. (b) Asample of bulk sequencing data includes multiple cells from differentclones. Their DNA content are extracted and are mixed together. Theassociation of cells’ genotype are lost in this step. (c) Horizontal blacklines represent reads from the bulk sample that are aligned to the ref-erence genome. Stars on the reads represents a variant allele comparedwith the reference genome. The blue star is a low prevalence muta-tion with a weak signal. (d) Single cells are separated, individuallysequenced, and separately aligned to the reference genome. Data hasthe cells’ genotype information.1.2 The landscape of genetic variants in tumoursOver the evolutionary process of cancer, tumour genome may acquire various ge-netic variants including single nucleotide variants (SNV) and structural variants(variants that are across a larger number of nucleotide sites).1.2.1 Single nucleotide variantsSNVs are point-like mutations that are a result of a single nucleotide substitution inthe cells genome. Somatic SNVs refer to variants that are present in tumour cells.Germline SNVs refer to variants that are present in both normal and tumour cells.Despite being a relatively simple genomic aberration, detection of somatic SNVs ina small sub-population of cells (rare clones) is challenging. For example, Figure1.1 panel c illustrates a high prevalence SNV and a low prevalence one by black andblue stars, respectively. It is obvious that while detection of high prevalence SNVsis rather trivial, detection of low prevalence ones is challenging due to weak signalto noise ratio. The weak signal is mainly because of (i) the presence of normalcells in the sample, (ii) copy number changes, and (iii) existence of SNVs only in asmall population of cells.51.2.2 Structural variantsIn contrast with SNVs that are point-like, structural variants (SV) are changes in alarge region of the genome. SVs are classified into different types depending on thenature of the variants. Deletion or multiplication of a region of the genome is calledcopy number variant (CNV). New genome arrangements resulting from genomebreakage at different positions and subsequent attachment of the broken ends isreferred as structural rearrangement. For example, inversion type rearrangementis when a segment of the genome is reversed end-to-end. Or translocation typerearrangement is when at least a segment of the genome moves to a new positionon the same or on another chromosome.1.2.3 SNVs and SVs in the study of tumour dynamicsBoth SNVs and SVs contain valuable information for the study of the evolutionarydynamics of tumour clones. The selective advantage and the abundance of SNVsacross the genome (numbering in the thousands) form a statistically robust marker.The signal of this marker is affected by CNVs. For example, the autosomal chro-mosomes of a diploid cell have two copies of each genome segment. If a SNVexists only in one of the copies, the number of variant nucleotides (mutated alleles)is half of the total number of alleles. However, if there is a CNV overlapping withthe SNV locus, the variant allele ratio at that locus can be different. Therefore, con-sideration of CNV data can be instrumental in detection of SNVs, particularly thelow prevalence ones. In both of the suggested models in this thesis (Chapter 2 andChapter 3) CNV information is incorporated to improve the accuracy of detectingSNVs.1.3 Next generation sequencing technologiesNext generation sequencing (NGS) significantly decreases the sequencing cost throughmassive parallel sequencing of millions of DNA fragments. Whole human genomewas first sequenced using NGS [25, 78].There are different NGS platforms, however their workflow can be summarizedas follows (see Figure 1.2):Library Preparation - DNA molecules of cells from a dissected sample are col-6Genomic DNAFragmented DNA Adaptors Sequencing library(a)(b) Flow cell Sequencing libraryClustersBridge amplification cycleACGTSequencing cycleDigital imageClusters1234Image processingCluster 1 > Read 1 > ACGGTACCluster 2 > Read 2 > CTGTGCCACluster 3 > Read 3 > GAGGTAGCluster 4 > Read 4 > TCCCTAGSequencing output fileACGGTCCTCCCTAGCGGTGCTCCCACGGTGCTCCCTACGCATTCGGTCCTCCCTAG---------------------------------ACGGTGCTCCCTAGGCATReadsReference genomeAlignment(c)Figure 1.2: (Caption next page.)7Figure 1.2: Workflow of next generation sequencing technology (a) NGS li-brary is prepared by fragmenting DNA molecules and ligating special-ized adaptors to both fragment ends. (b) The library is loaded into flowcell, and the fragmented DNAs are attached to the surface. Each frag-ment is amplified into a clonal cluster through a polymerase chain re-action. (c) The flow cell is imaged, and the emission from each clusteris recorded. The emission wavelength and intensity are used to iden-tify the bases. This cycle is repeated n times to create a read of lengthn. Reads are aligned to a reference genome using bioinformatics soft-ware. The aligned reads are used as raw data to identify genetic variants(modified figure from [59]).lected. DNA molecules are fragmented, tagmented and ligated by adaptor se-quences. In this step, the association of fragments and cells is lost.Cluster Generation - The library is loaded into a flow cell where fragmentedDNAs are amplified in a process called polymerase chain reaction (PCR). In thisstep each fragment is amplified into distinct clusters through a bridge amplificationprocess (see Figure 1.2 panel a).Sequencing - The flow cell is imaged and the emission from each cluster is recorded.The emission wavelength and intensity are used to identify the base. This cycle isrepeated n times to create a read of length n. This image-processing step is prone toerrors. Different parameters characterize the accuracy of this step. These parame-ters are utilized as a measure of the errors in the subsequent analysis (see Figure 1.2panel b).Alignment and Data Analysis - The reads are aligned to a reference genome us-ing bioinformatics software. Different post processing steps may be required forcorrecting sequencing and alignment errors. After alignment, the genetic variantscan be inferred from the reads by comparing them with a reference genome (seeFigure 1.2 panel c). The changes in genetic variants across time or across differentanatomical locations demonstrate the clonal dynamics of the tumour.In NGS data, each locus is covered by multiple overlapping reads that eachrepresents a DNA fragment from the original sample [6]. This data format providesa major benefit for the analysis of variants in the genome. The reason is that thedigital counts of the nucleotides (Adenine, Cytosine, Guanine, or Thymine) in each8locus quantify the nucleotides statistics. In the context of cancer, the allele countsare proportional to the population prevalence of an allele. This fact can be used toinfer the prevalence of the clonal populations [100].Several sequencing technologies were developed to study cancer genomes. Wereview two categories that are used to generate data for the models described inthis thesis.1.3.1 Bulk genome sequencingBulk sequencing technologies sequence multiple cells of a sample tissue altogether.Therefore, the output data is representative of the genome of a pool of tissue cells.The cells most likely belong to different populations (see Figure 1.1 panel a and b).The data must be deconvolved in order to identify the underlying genomic profilesof each cell or clone.The sequencing can be done across whole genome (whole genome sequencing)or only across the exome regions (exome sequencing). The sequencing can also bedone across a pre-specified segments of the genome (targeted sequencing).Whole genome and exome sequencingBoth whole genome and whole exome sequencing start with a large pool of cells.The cells are lysed, and their DNA contents are extracted. The cell genotype asso-ciation is lost. In whole genome sequencing (WGS) the entire genome of the ad-mixed DNA is sequenced, while in whole exome sequencing (WES) only the exomesegments of the genome are amplified and sequenced. In present whole genomesequencing, samples from DNA fragments are sequenced with the median coverageof 106x [32, 94, 96, 107, 108]. It means the average number of reads covering aposition of genome after alignment is about 106. The median coverage is expectedto increase with advances in the sequencing technologies. This number is largerin WES since it does not include all segments of the genome. One issue with bothWGS and WES is that the input DNA comes from a pool of cells. It means normalcells are sequenced together with cancerous (tumour) cells. Therefore, the resultantreads can represent DNA fragments from both the tumour cells and uninterestingnormal cells. In addition, the population of tumour cells is heterogeneous. All such9factors complicate detecting genetic variants (SNV or CNV) unique to rare clones.Stochastic sampling of DNA fragments and low coverage of WGS or WES makeDNA fragments of rare clones not represented or weakly represented by the reads(see the blue star mutation in Figure 1.1 panel (c). Copy number variants canalso affect the observed mutation signal. For example, in Figure 1.1 panel a, thegreen and purple mutations have similar prevalence. However, due to copy numbervariants the allelic prevalence of purple mutation is lower than the prevalence of thegreen one (see Figure 1.1 panel c). The CNV can make the SNV signal weak, andthat causes SNV detection tool to misinterpret it as an artifact. Therefore in bulksequencing data, the allelic prevalence does not directly and properly represent theclone populations. The issue of the weak signal of SNVs can be partly addressedby increasing the sequencing coverage rate.Deep-targeted sequencingTo overcome the issue with low coverage in WGS and WES, deep-targeted sequenc-ing can be used [32, 107, 108]. In deep-targeted sequencing only selected regionsof interest from the genome are amplified. The regions are usually the aberrantpositions that are initially identified from WGS. Their polymerase chain reaction(PCR) primers for the sequencing are then designed accordingly. The advantageof this approach is that only a small portion of the genome is sequenced. This re-sults in a data with a higher coverage of about 1000x-10000x. The higher coverageis beneficial through providing a stronger signal for accurate detection of SNVs(particularly the ones that are only present in a small sub-population of cells).Accurate detection of SNVs is critical for inferring the prevalence of underlyingsub-populations and their dynamics.1.3.2 Single cell sequencingSingle Cell Sequencing (SCS) technology sequences the genome of individual cells.This addresses the problem of DNA fragments being disassociated from the cells.In SCS, the first step is to isolate a cell from the primary sample to produce a viablesingle cell without a bias to any specific sub-population. The process is based onrecent developments in flow sorting or microfluidics [10, 89, 128, 141]. A normal10single cell contains only about 7 pg of DNA which is not enough material for cur-rent sequencing platforms. So, the cell needs to go through either whole genomeamplification (WGA), targeted polymerase chain reaction based protocols, or tag-mentation approaches to produce enough material for sequencing [76, 82, 141].SCS technology is expensive; it is in its early stages; and it is prone to differenttechnical errors. The errors include allele drop-out, uneven coverage across thegenome, and sequencing multiple cells (if two cells are sequenced together, themeasurement will be doublet). Figure 1.1 panel d, synthetically depicts the readsfrom single cell sequencing aligned to the reference genome. As shown, the datais sparse with uneven coverage across the genome.1.4 Literature review on SNV calling modelsAccurate detection of somatic single nucleotide variants is important in definingclonal composition of the human cancers. The genomic data from tumour tissueis accessible either through bulk sequencing of multiple cells or single cell se-quencing. However, each data type has its own challenges. The data from bulksequencing is generated from admixed DNA material of a large pool of cells. Thepooled data complicates the detection of SNVs that are only present in a subsetof cells (small clone). In single cell sequencing, the cell genotype association ismaintained. However, detection of SNVs remains challenging because the data isvery sparse with uneven coverage across the genome.1.4.1 Literature review on SNV calling from bulk sequencing dataTheoretically, SNVs can be identified given data with enough read depth coverageregardless of their variant allelic ratio. However, detecting SNVs with high confi-dence is challenging due to the artifacts and the limited read depth coverage. Mostof the artifacts are present in low frequency and they interfere with the detection oflow prevalence SNVs [135]. SNV calling has been an active research field for years.The state-of-the-art tools are designed to detect SNVs from matched tumour-normalsamples, a single tumour sample, or multiple normal and tumour samples. Matchednormal and tumour samples are required to detect somatic SNVs. If the matchednormal sample is not available or calling germline is desired, a single sample SNV11caller will be needed. In addition, there has been an increase in deciphering clonalpopulation dynamics over clonally related multiple samples (same patient samplestaken at different points in time or space). This has resulted in abundance of WGSor WES data from multiple samples. With the access to this data, methods incor-porating multiple clonally related samples are proposed to detect SNVs with higheraccuracy.Matched tumour-normal SNV callingThere is a wide range of methods proposed for SNV calling using matched normaltumour samples. Most of the approaches can be categorized into four groups: (i)Bayesian approaches, (ii) statistical test approaches, (iii) heuristic approaches, and(iv) machine learning approaches.Bayesian approach SomaticSniper, JointSNVMix2, Samtools, Virmid, FaSD-somatic and SNVSniffer are examples of the models taking Bayesian approach fordetecting SNVs [24, 67, 79, 81, 99, 127]. These methods evaluate the posteriorprobability of the joint genotype of normal and tumour samples (assuming theyare diploid). The methods may apply some filtering as the post processing stepto improve the performance. The posterior probability of the joint genotype iscalculated byP(GT ,GN |DT ,DN) = P(DT ,DN |GT ,GN)P(GT ,GN)∑gT ,gN∈G P(DT ,DN |gT ,gN)P(gT ,gN). (1.1)GT and GN , tumour and normal genotypes, take value from G = {AA, AC, AG,AT , CC, CG, CT , GG, GT , T T}, where A, C, G, T represents DNA nucleotides.DT and DN are tumour and normal data (aligned reads), respectively. The priorgenotype probability, P(GT ,GN) depends on the prior knowledge of mutation rate.The joint genotype likelihood follows a Binomial probability distribution with themean parameter that depends on the genotype and sequencing error rate [135].SomaticSniper and FaSD-somatic detect a site as a SNV if the probability of hav-ing different normal and tumour genotypes is high [79, 127]. SAMtools followsthe same idea but uses log-likelihood ratio instead of posterior probability [24].JointSNVMix2, Virmid and SNVsniffer simplify the 10 genotype state of G into12AA, AB and BB [67, 81, 99]. A denotes the allele matching the reference and B isthe one that does not match. Therefore, the probability of marking a site as a SNVis calculated byP(Somatic) = P(AA,AB)+P(AA,BB), (1.2)where the normal genotype alleles match the reference (AA or homozygous ref-erence). The tumour genotype is heterozygous in the first term (AB), and a non-reference homozygous genotype in the second term (BB) [135]. Virmid consid-ers the tumour as a mixture of normal and tumour samples with somatic muta-tions [67]. The model infers the tumour content and joint genotype. JointSNVMix2applies a hierarchical Bayesian model to estimate the joint genotype probabili-ties [99]. SNVsniffer initially calls high confidence SNVs using a set of heuristicthresholds. Then, the resultant low confidence calls are further tested using thejoint genotype model [81].The assumption that the normal and tumour samples are diploid, is relaxed inStrelka, MuTect, deepSNV, and EBCall [22, 44, 105, 111]. These models usuallyincorporate the joint variant allele frequencies (VAF), denoted by ( fT , fN), insteadof the joint genotypes (GT ,GN). Strelka infers the VAFs of normal and tumour sam-ples. Assuming the normal sample genotype as a homozygous reference, the modelcalculates the probability of difference in VAFs of normal and tumour samples.MuTect considers two models (i) wildtype model, and (ii) mutation model. Thewildtype model assumes that the non-reference reads come from artifacts, whereasmutation model assumes that the variant alleles come from an unknown frequency.MuTect then computes which of the aforementioned models are more likely [22].In local realignment, the reads are assembled together, and the read supportsare then counted. Mutect2, Platypus, FreeBayes incorporate the local realignmentstep [3, 39, 98]. This strategy improves the performance particularly in regionsprone to errors. It also gives information about the coexistence of mutations.Statistical test approach In this approach the model’s null hypothesis is wild-type and the alternative hypothesis is variant [19, 44, 111, 131]. LoLopicker andLoFreq use a sequencing quality score to infer the error rates [19, 131]. DeepSNVestimates the error rates directly from the input normal and tumour samples, whileEBcall uses other independent control samples [44, 111]. Note LoLoPicker and13deepSNV are particularly designed for deep targeted sequencing data to detect lowprevalence SNVs [19, 44].Heuristic approach The methods with heuristic approach make use of differentthresholds to identify mutation candidates, and then apply ad hoc rules to deter-mine somatic SNVs [53, 69, 74, 135]. For instance, VarScan2 requires the variantallelic ratio of more than 8% with at least 2 reads supporting the variant [69]. Thethresholds are mainly set above the expected artifact rate. For the loci passed thethresholds, VarDict, VarScan2 and Shimmer apply Fisher’s exact test in a two-by-two table of read counts [53, 69, 74]. Small p-value denotes somatic SNVs thathave normal and tumour variant reads significantly different.Machine learning approach Another group of methods for calling SNVs arebased on machine learning. The methods classify mutations into either (i) wild-type/germline or (ii) somatic. MutationSeq and SNooper use supervised classifierslike random forest, while BAYSIC applies an unsupervised latent class model todetect somatic SNVs [17, 31, 119].Single sample SNV callingSingle sample SNV calling is required when the matched normal sample is notavailable (e.g Formalin-Fixed Paraffin-Embedded (FFPE) samples), or when SNVcalling is required for a tumour deep targeted sequencing data. Outlyzer, ISOWN,and SomVarIUS, Pisces, SPLINTER, SNVer, SiNVICT are examples of the modelsparticularly designed for deep targeted sequencing data [35, 64, 70, 86, 116, 124,129]. Outlyzer, ISOWN, and SomVarIUS use FFPE data in their analysis [64, 86,116], and SNVMix2 use single tumour WGS/WES data [47].OutLyzer estimates the background noise and calls the variants that are abovethe noise level [86]. SPLINTER generates an error model of the pooled samplesto detect low prevalence variants. SNVer calls variants if the VAF exceeds a set ofthresholds based on Binomial distribution [129]. SNVMix2 model is based on theinference of the posterior probability of different genotypes [47].ISOWN, SiNVICT and SomVArIUS are other single sample SNV calling meth-ods that are capable of distinguishing somatic and germline SNVs [64, 70, 116].ISOWN uses a supervised machine learning approach to classify the variants. The14classification uses different features including the membership of databases (likeCOSMIC for somatic mutations and dpSNP for germlines) [64]. SomVarIUS con-siders the VAF distribution of heterozygous germline SNPs as a measure to detectSNVs. The loci that are placed on the left tail of the distribution, are identified assomatic SNVs [116]. SiNVICT is designed for detecting low prevalence SNVs fromcirculating tumour DNA (ctDNA). The loci that have VAF significantly lower than0.5, are called somatic SNVs [70].Multiple sample SNV callingIntra-tumour heterogeneity complicates the detection of SNVs. In particular, theclone specific mutations that are only present in a small fraction of cells may bemissed. Increasing interest in deciphering clonal population dynamics from mul-tiple samples collected over time or over anatomical space, results in availabilityof multiple clonally related sequencing data. With the access to this data, methodslike MultiSNV, SNV-PPILP and MultiGems are proposed [63, 87, 126]. MultiSNVjointly analyses all available samples under a Bayesian framework to improve theperformance of calling SNVs that are shared [63]. SNV calls from GATK are re-fined and corrected by using phylogeny information across multiple samples [83].MultiGeMS calls SNVs using a statistical model selection procedure and accountsfor enzymatic substitution sequencing errors [87].1.4.2 Literature review on SNV detection from single cell sequencingdataSingle cell sequencing data brings an opportunity for detecting SNVs at each cell.However, the data is not perfect and introduces challenges including: low breadthand depth coverage, biased allelic counts, and possibility of sequencing multiplecells (most probably two cells). Monovar addressed the low and uneven data cover-age by pooling the sequencing data across cells without assuming any dependencyacross sites [140]. SCcaller identifies SNVs while accounting for local allelic am-plification biases. However, it cannot recover mutations from drop-out events orfrom loss of heterozygosity [33]. Conbase and LiRA utilize read phasing to correctfor errors and allelic drop-out [12, 54]. Conbase unlike LiRA performs joint vari-15ant calling across the population of cells. Although Conbase and LiRA use readphasing in order to decrease false discovery rate (FDR) and increase specificity,the analysis of bases is only possible for regions in proximity to single nucleotidepolymorphisms (SNPs). This requirement refrains to call SNVs across any basesthat are covered by reads. SCIΦ detects SNVs while inferring the phylogenies withmutations propagated along tree branches [115]. Although SCIΦ claims to detectmutation in single cells with very low or even no coverage, it requires at least twocells to show an alternative nucleotide count of at least three [115]. Genotyper is amethod to infer clonal genotypes. It is based on a statistical model coupled with amean-field variational inference approach [101].The advent of DLP protocol resulted in data that consists of thousands of cellswith very low coverage (about 0.07x). The existing methods are only scalable tohundreds of cells and require more coverage for calling SNVs. Therefore, a SNVcaller that considers the characteristics and scale of the data is needed.1.5 Motivations and research contributionWe discussed the evolutionary process of cancer where mutations show a recordof ancestral relationships. Therefore, the accuracy of mutation calling is criticalin inferring the underlying phylogeny and the clonal dynamics. We reviewed theexisting sequencing technologies and the methods for detecting SNVs using (i) bulksequencing data, and (ii) single cell sequencing data. Conventionally, mutations arethe input for inferring the underlying clones’ phylogeny. In this thesis, we showthat injecting underlying clonal structure into the SNV detection model improvesthe performance of detecting SNVs. Our method improves the accuracy of inferringthe phylogeny and detecting mutations through multiple iterations.The first contribution is to use the underlying clonal information (can be in-ferred from deep targeted sequencing data) in a model to accurately detect SNVs.This is done using WGS or WES data from multiple clonally related bulk tumoursamples. The samples can be from different anatomical locations or multiple sam-ples from different time points. We encode detecting SNVs in a generative prob-abilistic framework, and we call it MuClone [34]. It performs joint statistical in-ference on multiple observations (from multiple samples) of the allele counts of a16locus of interest. We evaluated the performance of MuClone on both synthetic andreal data in Chapter 2.The second contribution is to use the underlying phylogeny (can be inferredfrom CNV data) in a model to detect SNVs in every individual cell. The modelutilizes single cell sequencing data from one or multiple clonally related tissuesamples. We call the model CellMutScope, and it incorporates a Bayesian sta-tistical framework. It performs joint statistical inference on the allelic counts ofmultiple cells at a locus of interest, and outputs the posterior probability of being amutation at the loci for every individual cell. Also, we frame incorporation of SNVdata in the Corrupt model (Corrupt infers the underlying phylogeny of cells usingCNV data). The extended tree with both SNV and CNV discloses the underlyingtree phylogeny in more detail. The performance of CellMutScope is benchmarkedis Chapter 3, and its performance on real data is presented in Chapter 4. Chapter 5concludes by a summary of contributions and a discussion of future works.17Chapter 2Somatic mutation detection andclassification throughprobabilistic integration of clonalpopulation structure“A very great deal more truth can become known than can beproven.”- Richard FeynmanThis chapter introduces MuClone, a new model for detection of SNVs across mul-tiple clonally related whole genome sequencing data or whole exome sequencingdata. In addition to detection of SNVs, MuClone classifies SNVs into biologicallymeaningful groups to allow study of clonal dynamics. The problem and the ex-isting approaches are reviewed in Sections 2.1 and 2.2, respectively. Detection ofSNVs using WGS or WES of multiple tumour samples by MuClone is presented inSection 2.3. We demonstrated that incorporating clonal information into joint anal-ysis of multiple samples improves SNV detection. This has a particular advantagefor low prevalence SNVs. The performance of MuClone on both synthetic and realdata is reported in Section 2.4. Section 2.5 concludes and discusses the next steps.182.1 IntroductionGenomic accumulation of somatic SNVs can disrupt the regular activity of cells andcan result in cancer initiation and progression. Collectively, the complete repertoireof SNVs across a cancer genome (numbering in the thousands) form a statisticallyrobust marker for inferring clonal populations to study tumour evolution. As such,accurate detection of all somatic SNVs, including those with low prevalence, is vi-tal as they can define clones with phenotypic properties of interest. Mechanistic as-sociation of specific clones with properties such as treatment resistance, metastaticpotential, and fitness under therapeutic selective pressures remains a key objectiveof biomedical investigators studying tumour progression.Phylogenetic analysis can encode the evolutionary lineage of tumour cells acrosstime and anatomic space [46, 60, 71, 84, 90, 107, 136]. Sequencing of multiplesamples of a cancer to reconstruct evolutionary patterns and drug response profilesis increasingly common. For example, in rapid autopsy programs, at the time ofa patient’s death, tens to hundreds of metastatic samples are collected for futurestudy [113, 134]. Recent multi-sample sequencing studies in renal, lung, ovary,breast, colorectal and other cancers have revealed striking evolutionary and clin-ically important properties of cancers [36, 42, 60, 84]. However, the analyticalmethods to detect mutations from such experimental designs are still immature,and few studies have leveraged shared statistical strength across samples to detectmutations with greater sensitivity.In the limiting case, all cells likely harbour unique genomes; however due tothe nature of branched evolutionary processes, clones can be coarsely modelledas major clades in the cell lineage phylogeny of a cancer. These clades share themajority of mutations, and therefore define first approximations to the genotypesof clones. Clonal genotypes and their relative abundances in the cancer cell pop-ulation can be approximated by clustering mutations measured in bulk tissues andestimating their cellular prevalences (the variant tumour cell fraction) [100, 139].Phylogenetic algorithms mostly use mutations (represented as binary geneticmarkers), as inputs to infer the branched evolutionary lineages of tumour cells [30,95]. Thus, mutation detection accuracy will ultimately impact the performance ofphylogenetic inference algorithms.19Detection of low prevalence mutations is a major challenge due to typicallysmall signal to noise ratios, owing to: (i) contamination by normal cells; (ii)genome copy number alteration; and (iii) the presence of mutations in only a smallfraction of tumour cells (intra-tumour heterogeneity). In this work, we illustratethat knowledge of the clonal population structure improves detection of mutationsdefining low prevalence clonal genotypes.2.2 Literature reviewAlthough SNV calling algorithms are ubiquitous in the literature, it remains chal-lenging to detect low prevalence mutations. Algorithms have been developed forcalling mutations from a single sample [47, 69], paired (matched normal and tu-mour) samples [22, 31, 67, 99, 105], or multiple samples [63, 126]. We list severalpopular algorithms. Mutationseq uses a feature based classifier for calling mu-tations [31], where the features are constructed from matched paired normal andtumour samples. Strelka is a method for somatic SNV and small indel detectionfrom sequencing data of matched normal and tumour samples [105]. It is basedon a Bayesian approach that uses normal and tumour samples’ allele frequencieswith normal expected genotype structure. MuTect uses a Bayesian classifier thatemploys various filters to ensure high specificity to detect mutations from matchednormal and tumour samples [22]. FreeBayes uses short read alignments to call themost likely genotypes for the population at each position. It can be run in singlemode using only one tumour sample or in multiple mode utilizing multiple tumoursamples from the same patient [39]. FreeBayes can detect somatic mutations ifgermlines are manually removed. MultiSNV jointly analyses all available sam-ples under a Bayesian framework to improve the performance of calling sharedmutations [63]. SNV calls from GATK [83] are refined and corrected by usingphylogeny information across multiple samples [103, 126].2.3 MethodMuClone uses previously known cellular prevalence information to improve mu-tation detection and classification. For each sample, MuClone detects mutationsfrom joint analysis of multiple samples. We encode this process in a generative20probabilistic framework to perform joint statistical inference of multiple observa-tions (from multiple samples) of the variant allele counts of a mutation of interest.The inputs to the model are: the number of variant reads, and the depth for a setof sequenced loci from multiple samples derived from the same patient; a measureof allele specific copy number at each locus, in each sample, with tumour content;and the cellular prevalence and the abundance of underlying mutation clusters. Mu-Clone outputs (i) a probability for each locus, at each sample, of being a mutation,and (ii) its cluster.The probabilistic graphical model of MuClone is depicted in Figure 2.1.2.3.1 Problem formulationIn our proposal, MuClone, we exploit prior knowledge of tumour content, tumourcellular prevalence, and copy number information across multiple samples to im-prove detection of somatic single nucleotide variants, and in particular, low preva-lence ones. Our model uses mutation clusters and copy number information ob-tained from standard approaches [50, 125]. In the first step, a set of stringent SNVcalls or validated SNVs (using targeted sequencing data) is used to infer muta-tion cluster information. Then, MuClone uses the inferred cluster information tomore accurately call mutations across genome (whole genome or exome sequenc-ing data). In addition to calling mutations, MuClone also classifies mutations intoclusters based on cellular prevalence. This provides the user with the opportunityto profile mutation changes across time and space, and adds a rich layer of inter-pretation into the detection process.2.3.2 Model descriptionWe first define gm,n, the genotype of a given locus n in sample m. Samples areindexed by 1 . . .M and loci are indexed by 1 . . .N. Given the copy number cm,n,the possible genotype states are G = {A . . .A,AA . . .B,A . . .BB, . . . ,B . . .B}, whereeach element has a length equal to cm,n. For example, the genotype ABB refers tothe genotype with one reference allele A and two variant alleles B. For simplic-ity, we assume the number of reads containing the variant alleles, bm,n at a givenlocus with genotype gm,n and read depth dm,n follows a Binomial distribution with21dm,n bm,n ψm,n pim,ntmsseqZnΩm = 1 . . .Mn = 1 . . .N1Figure 2.1: Probabilistic graphical model for MuClone: white nodes are un-observed variables; grey shaded nodes are observed variables. The vari-ables m ∈ {1 . . .M} and n ∈ {1 . . .N} index the samples and the loci,respectively. In sample m, dm,n is the total number of reads aligned atlocus n and bm,n is the number of aligned reads with B alleles. Thegenotype state is ψm,n and pim,n is the prior over the genotype states.The tumour content of sample m is tm and the error rate is εseq. The pa-rameter s stands for the precision parameter. Tumour clusters prior andtheir cellular prevalence information are encoded in Ω and the variableZn denotes the mutation cluster.genotype specific variant probability p(gm,n)bm,n|dm,n, p(gm,n)∼ Binomial(dm,n, p(gm,n)) . (2.1)22For gm,n ∈ G , the variant probability p(gm,n) : G → [0,1] is defined asp(gm,n) =B(gm,n)C (gm,n)B(gm,n) 6= 0, ,B(gm,n) 6= C (gm,n)εseq B(gm,n) = 0,1− εseq B(gm,n) 6= 0,B(gm,n) = C (gm,n),(2.2)where B(gm,n) : G → N and C (gm,n) : G → N return the number of the variantalleles of genotype gm,n and its copy number, respectively. For exampleB(ABB) =2 and C (ABB) = 3. The variable εseq > 0 is a small positive constant that accountsfor sequencing error. It allows for non-zero variant reads, due to sequencing error,when there are no variant alleles in genotype gm,n.However, since the sequenced reads are independently sampled from an infinitepool of DNA fragments, at a given locus, each read may belong to the normal,reference, or variant population. The normal population stands for normal cells; thereference population are tumour cells which do not have the mutation at the givenlocus; and the variant population are the ones carrying the mutation. Therefore,using a single genotype state, gm,n, introduces error into our analysis. To accountfor this fact, we consider using the full genotype state, ψm,n = (gNm,n,gRm,n,gV m,n),at a given locus n and sample m, to model the number of variant reads. Normalpopulation fraction is 1− tm where tm is the tumour content (the proportion oftumour cells in the sample) of sample m and the cellular prevalence of the mutationis φ zm which is the fraction of tumour cells carrying the mutation. According to ourprior knowledge, we assume mutations are clustered into K clusters. For a givenlocus n, Zn = z ∈ {1, . . . ,K} defines which cluster the mutation belongs to. If aposition is not a mutation then it belongs to wildtype cluster identified by Zn = 0.Therefore, for a mutation with cellular prevalence φ zm and tumour content tm,the variant allele probability is denoted by ξ (ψm,n,φ zm, tm). It is proportional to thesum of the (properly scaled) variant probabilities from each population:ξ (ψm,n,φ zm, tm) ∝ (1− tm)C (gNm,n)p(gNm,n)+tm(1−φ zm)C (gRm,n)p(gRm,n)+tmφ zmC (gV m,n)p(gV m,n), (2.3)23where the first term (1− tm)C (gNm,n)p(gNm,n) is proportional to the probabilityof sampling a read containing variant allele from the normal population, and thesecond and third terms, tm(1−φ zm)C (gRm,n)p(gRm,n) and tmφ zmC (gV m,n)p(gV m,n),are proportional to the probabilities of sampling a read containing variant allelesfrom the reference and variant populations, respectively.Considering the full genotype state, the number of reads containing the variantalleles at a given locus n that belongs to cluster Zn follows a Binomial distributionwith probabilityµ(Zn) =εseq if Zn = 0ξ (ψm,n,φ zm, tm) if Zn = z and z ∈ {1, . . . ,K}, (2.4)where εseq accounts for sequencing error in wildtype cluster and ξ (ψm,n,φ zm, tm)is the variant alleles probability for nth locus, mth sample from zth cluster. Ac-cording to Equation (2.3), tumour content and cellular prevalence information areincorporated to estimate ξ (ψm,n,φ zm, tm).Since empirical evidence shows that variant read data is overdispersed, we re-place the Binomial model (2.1) with a BetaBinomial modelbm,n|dm,n,µ(Zn),s∼ BetaBinomial(bm,n|dm,n,µ(Zn),s) , (2.5)where µ(Zn) is the expected variant alleles probability and the hyperparameter sis the precision parameter of the BetaBinomial distribution. The BetaBinomialdistribution in Equation (2.5) assigns a small chance for mutation when the locusis wildtype, otherwise it is governed by the prior clonal information.To fully express our model, for each locus, we assume the genotype state fol-lows a categorical distribution with probability vector pim,n ∈ [0,1]|G | whose ithelement is the probability of the ith genotype state,ψm,n|pim,n ∼ Categorical(pim,n) . (2.6)The number of possible genotype states, denoted by |G |, is finite given the copynumber information. For simplicity, we assume every element of pim,n is equal to1|G | .24In addition, we also assume that the clonal assignment of a locus, denoted byZn, follows categorical distribution with probability vector τ :Zn|τ ∼ Categorical(τ ) . (2.7)Our probabilistic framework can be succinctly written asbm,n|dm,n,µ(Zn),s∼ BetaBinomial(bm,n|dm,n,µ(Zn),s) ,ψm,n|pim,n ∼ Categorical(pim,n) ,Zn|τ ∼ Categorical(τ ) .(2.8)2.3.3 InferenceBased on the generative model introduced in (2.8) mutations are inferred via theposterior probability distribution of a locus n belonging to cluster z:P(Zn = z|bm,n,dm,n,s) ∝ τzM∏m=1∑i∈Ipim,niL (Zn = z|bm,n,dm,n,s), (2.9)where the variable i indexes pim,n over the genotype states, I = {1 . . . |G |}. The pos-terior probability of locus n belongs to cluster z is proportional to the likelihood ofobserving bm,n number of nucleotides matching the variant alleles times the priorover tumour cluster z. The tumour cluster prior τz is the fraction of mutations be-longing to cluster z and it has been tuned according to wildtype prior; the tumourcluster prior and the cellular prevalence information are encoded in Ω. Wildtypeprior is our prior information if a locus is a mutation. If we don’t have any infor-mation we can set it to 0.5. The likelihood function,L (Zn = z|bm,n,dm,n,s), is theBetaBinomial distribution defined in 2.5.Based on basic decision theory, a decision can be extracted from a posteriordistribution given a loss function. Under the loss function `(z,z′) = 1[1[z = 0] 6=1[z′ = 0]], the decision is simply the maximum a posteriori (MAP). That is, if theprobability η of belonging to any of the tumour clusters is greater than 0.5, weconclude that the locus is mutated in at least one of the M samples. The value of η25isη =K∑z=1P(Zn = z|bm,n,dm,n,s).If locus n is mutated in at least one of the M samples, then the probability ofmutation, in each sample, is calculated separately asPm,n(mutant) = ∑j∈J∗mP(Zn = j|bm,n,dm,n,s), (2.10)where J∗m is the set of clusters of sample m whose cellular prevalences are greaterthan a fixed positive threshold called ΦT ,J∗m = {j | φ jm >ΦT}. (2.11)The threshold ΦT distinguishes the clusters of sample m in which their non-zerocellular prevalence are due to actual variant alleles. The default value of ΦT iszero. However, depending on the method used for estimating cellular prevalences,it can be set to another positive value, if some non-zero input cellular prevalencesindicate wildtype clusters.In addition, MuClone assigns the locus to cluster z∗ that maximizesz∗ = argmaxz∈{1,...,K}P(Zn = z|bm,n,dm,n,s). (2.12)This classifies mutations to one of the previously known clusters. The classificationof mutations helps in biological interpretation and phylogenetic analysis of thedata.2.4 Experimental results2.4.1 Synthetic dataIn this section, we examine the performance of MuClone on simulated data. Inwhat follows, we generate N loci from M samples with K underlying tumour mu-tation clusters with sequencing error rate εseq and tumour content tm.We first randomly generate an evolutionary relationship between clusters, viewed26S3 S2 S1 S00.00.20.40.60.81.0Prevalence 012345678910Cluster’s nameSamplesFigure 2.2: Distribution of the cellular prevalence of clusters across differentsamples in multiple runs.as a binary phylogenetic tree. Each node in the tree represents a mutation cluster.The root node represents the ancestral cluster. The cellular prevalences of the firstdescendant, φ1st , is sampled from a Uniform distribution over [0,φparent], whereφparent is the cellular prevalence of the parent node (cluster). The cellular preva-lences of the second descendant, φ2nd , is sampled from a Uniform distribution over[0,φparent−φ1st ], defined so that the sum of the children’s prevalences do not exceedtheir parent’s cellular prevalence. The absence or presence of each cluster, in eachsample, is sampled from a Bernoulli distribution that assigns equal probability toboth outcomes. Distribution of the cellular prevalence of 10 clusters across 4 sam-ples in multiple runs is depicted in Figure 2.2. Figure 2.3 is an example of cellularprevalences of clusters across different samples in one random run.If a cluster is not present in a sample, the corresponding cellular prevalencewill be 0. See Figure 2.4 for an example of this process. Loci are assigned to a27S3 S2 S1 S0Samples0.00.20.40.60.81.0PrevalenceyCluster name012345678910Figure 2.3: Distribution of the cellular prevalence of clusters across differentsamples in one random run.cluster uniformly at random from {0, . . . ,K}, where cluster 0 represents the wild-type cluster and {1, . . . ,K} are mutation clusters. For each locus, in each sample,the number of reads overlapping the locus (depth) is sampled from a Poisson dis-tribution with mean dm. Wildtype copy number is deterministically set to 2, and acopy number profile (major and minor copy number) is generated according to thefollowing steps: The total copy number, ct , is sampled uniformly at random from{1, . . . ,cmax}. An integer number, cb, is randomly (following a discrete Uniformdistribution) picked from 1 to ct , and ca is defined as ca = ct−cb. Lastly, the majorcopy number is set to the maximum of cb and ca; the minor copy number is set tothe minimum of those two values. Then, corresponding to each cluster, the num-ber of variant reads are sampled from the Beta-Binomial distribution described in28(b)(a)Mutation clusters treeS0 S1 S2 S3Sample0.00.20.40.60.81.0Prevalence 1, n=20%2, n=22%3, n = 30%4, n=2%5, n=25%Figure 2.4: A clonal information example. (a) Cellular prevalences of under-lying mutation clusters across multiple samples. (b) The relationshipbetween mutation clusters is represented as a tree.Equation (2.5) with precision parameter equal to 1000.Synthetic data evaluation We simulated 10 synthetic data for 20000 loci from4 samples of a patient, with 5 underlying clusters, including an ancestral cluster.The maximum copy number was 5, and error rate was 0.01. The average sequenc-ing depth was assumed to be 100 for all samples.To assess the performance and robustness of MuClone, we systematically shieldMuClone from clonal information (Figure 2.5). In particular, the cellular preva-lence information was perturbed by (i) adding noise to its value, or (ii) removingthe cellular prevalence information of the clusters. The noise was generated froma normal distribution with mean zero and standard deviations: 0, 0.01, 0.1, and0.2. The noise value, ν , was added to the cellular prevalence of the cluster, whilebounding the resulting value between 0 and 1, that is,φ ∗zm = min(max(φzm+ν ,0),1),where φ ∗zm and φ zm are the perturbed and original cellular prevalence of cluster z andsample m, respectively. The clusters which their clonal information was removed,were randomly chosen with equal probabilities. As expected, both sensitivity and29Figure 2.5: MuClone’s performance with inaccurate clonal information: 10synthetic datasets generated for 20000 loci, from 4 samples of a hypo-thetical patient, with 5 underlying clusters. The maximum copy numberis 5, error rate is set to 0.01, and average sequencing depth is approx-imately 100. To assess performance, we add noise from a mean zeronormal with standard deviation equal to 0 (dark purple), 0.01 (light pur-ple), 0.1 (light blue), and 0.2 (dark blue) to the cellular prevalence andalso remove the clonal information of different number of clusters. (a)Sensitivity and (b) Specificity of MuClone with parameters: wildtypeprior = 0.5, ΦT = 0.02, error rate = 0.01, tumour content = 0.75, andprecision parameter = 1000.specificity were highest with complete and accurate clonal information; see Fig-ure 2.5. This suggests that incorporating clonal information can improve mutationdetection accuracy and gives evidence to support MuClone’s approach. Further-more, since the sensitivities were only marginally impacted by adding noise to theclonal information, MuClone should be able to cope with modest misspecificationof the prior. However, specificity can decrease if the cellular prevalence is reducedto levels associated with the wildtype cluster and sensitivity can improve if addingnoise increases the cellular prevalence to levels associated with a removed mutationcluster.Naturally, accuracy was most severely impacted with reduced/corrupted clonalinformation; see Figure 2.5. For modest level of noise (noise standard deviation 030and 0.01), the sensitivity and specificity of removing various numbers of clusterswere compared through a Kruskal-Wallis test ( 4e−5 ≤ p-values ≤ 1e−4) whichshows that the change in performance due to clonal information is significant. Innoiseless settings, the confidence interval for the difference (of zero and four re-moved clusters) in mean sensitivity and specificity are [0.16,0.26] and [0.08,0.32],respectively. When the noise standard deviation is equal to 0.01, these intervals are[0.11,0.21] and [0.09,0.36].We also explored how sensitivity and specificity changes as a function of thewildtype prior and the threshold ΦT used to distinguish the cellular prevalence cut-off of a mutation cluster. In Figure 2.6, we tested MuClone with wildtype priorvalues 0.5, 0.75, and 0.99, and with ΦT values 0.001, 0.01, 0.02, 0.03, and 0.04.In the case that the wildtype prior equals 0.5, we assumed that a locus is equallylikely to be a mutation or not (when no other information is provided). MuClone’ssensitivity and specificity were near 1 for ΦT = 0.02 and wildtype prior equal to0.5. As expected, with small values of ΦT , the sensitivity and specificity decreasedsince it is difficult to distinguish between wildtypes and mutations with small cel-lular prevalences. The sensitivity also decreased for large values of ΦT becausemutations were miscalled as wildtypes. When the error rate was 0.01, and wild-type prior was 0.5, the optimal ΦT was about 0.02. We used these values for thefollowing experiments.The performance of MuClone was tested with various tumour content (from0.1 to 0.99) and different error rates (0.01 and 0.001); see Figure 2.7. For sam-ples with tumour content greater than 0.5, the sensitivity and specificity remainedclose to 1. Sensitivity and specificity decreased to only about 0.9 when the tumourcontent in the sample was as low as 0.1. These results establish promising per-formance over different ranges of tumour content with different error rates (likelyscenarios in real data). In addition, we also explored the performance of MuClonefor samples with different coverage (mean depth): 30, 60 and 100; see Figure 2.8.Intuitively, the performance was higher when we had more coverage. Since Mu-Clone leverages cellular prevalence information to improve the performance ofmutation detection, the performance gain was noticeable when the variant allelicratio resolution supports the given cellular prevalence resolution (in our analysisthe cellular prevalence of mutations were greater than 0.02).310.5 0.75 0.99Wildtype prior0.700.750.800.850.900.951.00Sensitivity0.5 0.75 0.99Wildtype prior0.700.750.800.850.900.951.00Specificity(a) (b)Figure 2.6: MuClone’s performance as a function of wildtype prior: 10 syn-thetic datasets generated for 20000 loci, from 4 samples of a hypothet-ical patient, with 5 underlying clusters. The maximum copy number is5 and error rate is set to 0.01, and average sequencing depth is approx-imately 100. We asses the performance for ΦT equal to 0.001 (darkpurple), 0.01 (light purple), 0.02 (white smoke), 0.03 (light blue), 0.04(dark blue). (a) Sensitivity and (b) Specificity of MuClone with param-eters: error rate = 0.01, tumour content = 0.75, and precision parameter= 1000.Figure 2.9 demonstrates how well mutations were classified by MuClone. Theinput clonal information was perturbed by adding noise from zero mean normaldistribution with with standard deviation 0.01 to simulate a more realistic scenario.In Figure 2.9(a), each bin (i, j) shows the fraction of mutations in cluster i that wereclassified into cluster j by MuClone. Figure 2.9(a) shows that 85% of mutationswere classified into the correct cluster.In order to show that the classification errors occurred between clusters withsmall phylogenetic distance, we define a misclassification index to quantify phylo-genetic distance; calculated asMisclassification index =∑i 6=j q(i,j)× dist(i,j)−distminidistmaxi −distmini∑i 6=j q(i,j), (2.13)320.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.99Tumour content0.700.750.800.850.900.951.00Sensitivity0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.99Tumour content0.700.750.800.850.900.951.00Specificity(a) (b)Figure 2.7: MuClone’s performance as a function of tumour content: 10 syn-thetic datasets generated for 20000 loci, from 4 samples of a hypothet-ical patient, with 5 underlying clusters. The maximum copy numberis 5, and average sequencing depth is approximately 100. We assesthe performance for error rate equal to 0.001 (light purple) and 0.01(light blue), and tumour content from 0.1 to 0.99. (a) Sensitivity and(b) Specificity of MuClone with parameters: wildtype prior = 0.5, φT= 0.02, and precision parameter = 1000.where q(i,j) is the number of mutations in cluster i that have been classified intocluster j, and the Euclidean distance between the cellular prevalences of clusteri and j is denoted by dist(i,j). The distance of the closest and farthest cluster tocluster i is denoted by distmini and distmaxi , respectively. In Figure 2.9(b), small mis-classification indices demonstrate that misclassified mutations occur between closeclusters. This can be interpreted as phylogenetically recently separated clusters.2.4.2 Real dataTwo real data sets with multiple samples for each patient were used to evaluate theperformance of MuClone. The first data set was multiple whole genome sequenc-ing data from 7 patients with high grade serous ovarian cancer. The second dataset was multiple whole exome sequencing data from 8 patients with non-small-celllung cancer (NSCLC).3330 60 100Depth0.8000.8250.8500.8750.9000.9250.9500.9751.000Sensitivity30 60 100Depth0.8000.8250.8500.8750.9000.9250.9500.9751.000Specificity(a) (b)Figure 2.8: MuClone’s performance as a function of different depth: 10 syn-thetic datasets generated for 20000 loci, from 4 samples of a hypotheti-cal patient, with 5 underlying clusters. The maximum copy number is 5,error rate is set to 0.01 and average sequencing depth is approximately100. (a) Sensitivity and (b) Specificity of MuClone with parameters:wildtype prior = 0.5, φT = 0.02, error rate = 0.01, tumour content =0.75, and precision parameter = 1000.High grade serous ovarian cancer We tested MuClone’s performance onwhole genome sequencing data (with depth 30x) from multiple tumour samplessurgically resected from high grade serous ovarian cancer patients [84]. The sam-ples were obtained from different spatially distributed metastatic sites. Brief de-tails about the number of samples for each patient, sample sites and the number ofvalidated loci for each patient are shown in Table 2.1. Germline mutations wereexcluded from the list.The copy number, tumour purity, and mutation cluster information for exper-imentally re-validated mutation status were taken from the phylogenetic study ofhigh-grade serous ovarian cancer (see the supplementary note of the paper [84]).Mutation clusters were estimated with PyClone [100] on the deep targeted se-quencing data (>1000x coverage) from the same samples and in three patientswith accompanying single cell sequencing data (see Table S16 in the phylogeneticstudy of high-grade serous ovarian cancer paper [84]). Copy number and tumour34-1 1 2 3 4 5True clusterWildtype clusterCluster oneCluster twoCluster threeCluster fourCluster fivePredicted cluster(a)1 2 3 4 5 6 7 8 9 10Different runs0.0000.0250.0500.0750.1000.1250.1500.1750.200Misclassification index(b)0.000.040.080.120.16Figure 2.9: MuClone’s classification performance: 10 synthetic datasets gen-erated for 20000 loci, from 4 samples of a hypothetical patient, with 5underlying clusters. The maximum copy number is 5, error rate is set to0.01. (a) Bin (i, j) shows the fraction of mutations in cluster i that wereclassified into cluster j by MuClone. 85% of the mutations are classi-fied correctly. (b) Misclassification index for 10 independent samples.MuClone with parameters: wildtype prior = 0.5, ΦT = 0.02, error rate= 0.01, tumour content = 0.75, and precision parameter = 1000.purity estimates were calculated with the TITAN software [50]. In order to elim-inate germlines, loci with variant nucleotides in the corresponding normal samplewere removed from the dataset. Then, the performance of MuClone was bench-marked against Strelka [105] (v2.0.15), MutationSeq [31] (v4.3.7), MuTect [3],FreeBayes [39] (v1.2.0-2), MultiSNV [63] and naive MuClone. Naive MuClone isa version of MuClone where no clonal information is provided (that is, all muta-tions are from an ancestral cluster).In Figure 2.10, the performance of MuClone is compared with other methodsexecuted with default settings. For each patient, p, we assessed performance byaveraging Youden’s index, sensitivity, and specificity across different samples. For35Patient Samples #Validated positions Anatomic samples1 6 153 Right Ovary Site 1-4; Omentum Site 1; Small Bowel Site l2 4 46 Omentum Site 1,2;Right Ovary Site 1,23 4 99 Right Ovary Site 1,2;Omentum Site 1; Left Ovary Site 24 5 69 Right Ovary Site 1-4;Right Pelvic Side Wall7 3 59 Left Ovary Site 1; Brain Metastasis; Right Pelvic Mass9 5 72 Right Ovary Site 1; Left Ovary Site 1; Omentum Site 1,210 4 136 Right Ovary Site 1-4Table 2.1: Summary of high grade serous ovarian cancer data set [84].patient p, with np samples, these are calculated asSensitivityp =1npnp∑i=1Sensitivityip,Specificityp =1npnp∑i=1Specificityip, (2.14)Youden’s indexp =1npnp∑i=1(Sensitivityip+Specificityip−1),where Sensitivityip, Specificityip and Youden’s indexip are the sensitivity, specificityand Youden’s index of sample i and patient p, respectively. In aggregate, Mu-Clone outperforms other methods by improving sensitivity without compromisingspecificity; see Figure 2.10. The performance (sensitivity, specificity and Youden’sindex) and the Receiver Operating Characteristic (ROC) curves for each patient areseparately depicted in Figures 2.11 to 2.24. False negatives arise mainly becausethe WGS data is under-represented (the average depth of the WGS data is about30x) and lacks variant alleles that are present in the targeted sequencing data. Falsepositives arise due to erroneous signal from sequencing technical artefacts.In Figure 2.10, Strelka, MutationSeq, MuTect and Naive MuClone have lowerperformance as they do not incorporate information across multiple samples. Free-Bayes was run on multiple samples and germlines were removed manually, butsince the method only considers tumour samples, it had the lower performanceversus other methods.To assess the performance of MuClone and MultiSNV, we conducted a two36MuCloneMultiSNVMutectNaive MuCloneStrelkaMutationSeqFreeBayes(a)0.20.40.60.81.0PerformanceYouden's IndexMuCloneMultiSNVMutectNaive MuCloneStrelkaMutationSeqFreeBayes(b)SensitivityMuCloneMultiSNVMutectNaive MuCloneStrelkaMutationSeqFreeBayes(c)SpecificityFigure 2.10: Performance comparison of different methods on whole genomesequencing data from patients with high grade serous ovarian cancer:(a) Youden’s index, (b) Sensitivity, and (c) Specificity across differ-ent mutation detection methods (from left to right: MuClone (darkblue), MultiSNV (orange), MuTect (light blue), Naive MuClone (yel-low), Strelka (purple), MutationSeq (brown), and FreeBayes (pink)).MuClone parameters are: wildtype prior = 0.5, ΦT = 0.02, tumourcontent = 0.75, error rate = 0.01, and precision parameter = 1000.sided t-test for the difference in the mean of Youden’s index evaluated on mutationcalls from MuClone and MultiSNV. The 95% confidence interval is [0.03,0.1],with p-value equal to 0.0006; this shows that the difference is statistically signif-icant. Importantly, MuClone improves sensitivity, enabling the detection of moremutations across the whole genome.Figure 2.25 depicts the classification of mutations into clusters relative to theground truth, as defined by running PyClone on the data (omitting singleton clus-ters [84]). Each bin (i, j) of Figure 2.25 shows the fraction of mutations in cluster ithat were classified into cluster j by MuClone; 93% are correctly classified by Mu-Clone. Moreover, we notice that misclassified mutations were classified into phy-logenetically similar clusters (the misclassification index for patient 1 was 0.015).Non-small-cell lung cancer We tested MuClone’s performance on early-stageNSCLC samples from the TRACERx data set [60]. (See Table 2.2). To help ob-tain the clonal and subclonal census, multiple tumour regions for each patient weresequenced by Illumina HiSeq. We used the copy number, purity estimate, andthe mutation cluster information available in the Supplementary Material of thepaper [60]. In the TRACERX study, the cellular prevalence was calculated from37MuCloneMultiSNVMutectNaive MuCloneStrelkaMutationSeqFreeBayes0.30.40.50.60.70.80.91.0PerformanceYouden’s IndexMuCloneMultiSNVMutectNaive MuCloneStrelkaMutationSeqFreeBayesSensitivityMuCloneMultiSNVMutectNaive MuCloneStrelkaMutationSeqFreeBayesSpecificity(a) (b) (c)Figure 2.11: Performance comparison of different methods on whole genomesequencing data for patient 1 with high grade serous ovarian cancer:(a) Youden’s index, (b) Sensitivity, and (c) Specificity across differ-ent mutation detection methods (from left to right: MuClone (darkblue), MultiSNV (orange), MuTect (light blue), Naive MuClone (yel-low), Strelka (purple), MutationSeq (brown), and FreeBayes (pink)).MuClone parameters are: wildtype prior = 0.5, ΦT = 0.02, tumourcontent = 0.75, error rate = 0.01, and precision parameter = 1000.the whole exome sequencing data on a set of stringent mutations that were se-lected from MuTect and VarScan2 results with post-processing. In addition, theTRACERx study added a few mutations to reduce missed subclonal mutations; seeSupplementary Appendix of TRACERx study [60].To compare the performance of MuClone with Strelka, MultiSNV, and MuTect,we randomly selected 8 patients with subclonal mutations from the TRACERx dataset (see Table S2 Supplementary Appendix 1 of the paper [60]). The TRACERxstudy generated a re-validated and curated list of mutations for their analysis; seeSupplementary Appendix 2 of TRACERx study [60]. The mutations with full copynumber information across all 8 patients were used as ground truth to evaluateperformance.We evaluated the false negative rates of mutation calling across several meth-ods; see Table 2.3. Altogether, out of 7238 mutations, MuClone missed 475 mu-tations while Strelka, MultiSNV and MuTect missed 7205, 5720, and 1086 mu-tations, respectively. Hence, borrowing statistical strength, as done in MuClone,across samples likely increases sensitivity to real mutations.380.0 0.5 1.0False positive rate0.00.51.0True positive rateROC curve Sample P1a, AUC=0.98139P1c, AUC=0.99925P1e, AUC=0.98952P1g, AUC=0.96773P1i, AUC=1.0P1k, AUC=0.99612Figure 2.12: MuClone’s Roc curves and the area under the curve (AUC) forpatient 1. MuClone parameters are: wildtype prior = 0.5, ΦT = 0.02,tumour content = 0.75, error rate = 0.01, and precision parameter =1000.We next ran MuClone, MultiSNV and MuTect on the whole exome data frommultiple samples of 8 patients to ascertain specificity. We note that MuClone re-moves reads with mapping quality less than 5 and for positions that have (i) avariant nucleotide in a normal sample, (ii) more than 40% filtered basecalls (Abasecall is filtered if more than 3 mismatches occur between the read and the ref-erence within a window of 20 bases on each side of the site.); or (iii) more than75% of the reads that cross the site have deletions in any of the samples [105].For exome sequencing data, mutations were called if the corresponding MuClonemutation probability is greater than 0.9. The other methods were executed with de-fault settings. The total number of calls and the number of common calls betweendifferent methods (restricted to positions with copy number information) at the pa-tient level is depicted in Figure 2.26. A high degree of variation across callers wasobserved.Altogether, MuClone called 13556 mutations while MultiSNV and MuTectcalled 31374 and 11915, respectively. MultiSNV output the largest number of39MuCloneMultiSNVMutectNaive MuCloneStrelkaMutationSeqFreeBayes0.30.40.50.60.70.80.91.0PerformanceYouden’s IndexMuCloneMultiSNVMutectNaive MuCloneStrelkaMutationSeqFreeBayesSensitivityMuCloneMultiSNVMutectNaive MuCloneStrelkaMutationSeqFreeBayesSpecificity(a) (b) (c)Figure 2.13: Performance comparison of different methods on whole genomesequencing data for patient 2 with high grade serous ovarian cancer:(a) Youden’s index, (b) Sensitivity, and (c) Specificity across differ-ent mutation detection methods (from left to right: MuClone (darkblue), MultiSNV (orange), MuTect (light blue), Naive MuClone (yel-low), Strelka (purple), MutationSeq (brown), and FreeBayes (pink)).MuClone parameters are: wildtype prior = 0.5, ΦT = 0.02, tumourcontent = 0.75, error rate = 0.01, and precision parameter = 1000.Patient Samples Cancer typeCRUK0003 4 AdenocarcinomaCRUK0004 4 AdenocarcinomaCRUK0005 4 AdenocarcinomaCRUK0013 5 AdenocarcinomaCRUK0062 7 Squamous-Cell CarcinomaCRUK0063 5 Squamous-Cell CarcinomaCRUK0065 6 Squamous-Cell CarcinomaCRUK0094 4 OtherTable 2.2: Summary of the NSCLC data set [60].400.0 0.5 1.0False positive rate0.00.51.0True positive rateROC curve Sample P2a, AUC=0.94827P2c, AUC=0.95039P2e, AUC=1.0P2g, AUC=1.0Figure 2.14: MuClone’s Roc curves and the area under the curve (AUC) forpatient 2. MuClone parameters are: wildtype prior = 0.5, ΦT = 0.02,tumour content = 0.75, error rate = 0.01, and precision parameter =1000.Patient MuClone MultiSNV MuTect StrelkaCRUK0003 52 350 60 430CRUK0004 16 188 36 240CRUK0005 212 1736 236 2040CRUK0013 26 490 270 540CRUK0062 30 469 42 609CRUK0063 75 445 40 510CRUK0065 60 1902 342 2640CRUK0094 4 140 60 196Table 2.3: Total number of false negative calls across multiple samples ofnon-small cell lung cancer patients for different algorithms.41MuCloneMultiSNVMutectNaive MuCloneStrelkaMutationSeqFreeBayes0.40.50.60.70.80.91.0PerformanceYouden’s IndexMuCloneMultiSNVMutectNaive MuCloneStrelkaMutationSeqFreeBayesSensitivityMuCloneMultiSNVMutectNaive MuCloneStrelkaMutationSeqFreeBayesSpecificity(a) (b) (c)Figure 2.15: Performance comparison of different methods on whole genomesequencing data for patient 3 with high grade serous ovarian cancer:(a) Youden’s index, (b) Sensitivity, and (c) Specificity across differ-ent mutation detection methods (from left to right: MuClone (darkblue), MultiSNV (orange), MuTect (light blue), Naive MuClone (yel-low), Strelka (purple), MutationSeq (brown), and FreeBayes (pink)).MuClone parameters are: wildtype prior = 0.5, ΦT = 0.02, tumourcontent = 0.75, error rate = 0.01, and precision parameter = 1000.calls in all of the samples, while MuTect and MuClone output similar number ofcalls. Figure 2.26 also demonstrates the mutations used in the TRACERx studyand their overlap with the mutation calls in different methods. The set of mutationsoverlapping between MuClone and TRACERx is most similar; this suggests thatthe increase in sensitivity conferred by MuClone does not come at the expense ofspecificity.We also explored the performance of MuClone when clonal information dif-fers in the number of input clusters or the value of the cellular prevalence; see Fig-ure 2.27. We perturbed the value of the cellular prevalences (estimated by PyClone)by adding noise from a mean zero normal distribution with different standard de-viations: 0, 0.01, 0.1, and 0.2. We see that MuClone is robust to slight changesof cellular prevalence values. We also shielded MuClone from different fractionsof the clonal information and that decreased performance more than adding noise.In general, this result shows that more accurate clonal information provides bettermutation detection.420.0 0.5 1.0False positive rate0.00.51.0True positive rateROC curve Sample P3a, AUC=0.99917P3c, AUC=0.97942P3e, AUC=0.97799P3g, AUC=0.98231Figure 2.16: MuClone’s Roc curves and the area under the curve (AUC) forpatient 3. MuClone parameters are: wildtype prior = 0.5, ΦT = 0.02,tumour content = 0.75, error rate = 0.01, and precision parameter =1000.2.5 ConclusionWe studied the use of clonal information for the purpose of somatic mutation de-tection and classification in multi-sample whole genome sequencing data. The pro-posed statistical framework uses the clusters cellular prevalences and copy numberinformation for detection and classification of low prevalence mutations. Our pro-posal, MuClone, outperformed other popular mutation detection tools, while pro-viding the added benefit of classifying whole genome sequencing mutations intobiologically relevant groups. Both synthetic and real data results showed that usingthe cellular prevalences of tumour clusters can improve mutation detection sensi-tivity. Importantly, our results suggest improvement in sensitivity can be achievedwithout compromising specificity.Since the accuracy of detecting mutations can affect the performance of phy-logenetic analysis, we suggest improvement from using MuClone will impact thefield of multi-region sequencing for cancer evolution studies. As the field matures,43MuCloneMultiSNVMutectNaive MuCloneStrelkaMutationSeqFreeBayes0.50.60.70.80.91.0PerformanceYouden’s IndexMuCloneMultiSNVMutectNaive MuCloneStrelkaMutationSeqFreeBayesSensitivityMuCloneMultiSNVMutectNaive MuCloneStrelkaMutationSeqFreeBayesSpecificity(a) (b) (c)Figure 2.17: Performance comparison of different methods on whole genomesequencing data for patient 4 with high grade serous ovarian cancer:(a) Youden’s index, (b) Sensitivity, and (c) Specificity across differ-ent mutation detection methods (from left to right: MuClone (darkblue), MultiSNV (orange), MuTect (light blue), Naive MuClone (yel-low), Strelka (purple), MutationSeq (brown), and FreeBayes (pink)).MuClone parameters are: wildtype prior = 0.5, ΦT = 0.02, tumourcontent = 0.75, error rate = 0.01, and precision parameter = 1000.we expect that the method presented here will be incorporated into more analyti-cally comprehensive modelling of whole genome sequencing data when multiplesamples are used to infer properties of clonal dynamics. Next steps are in devel-oping a unified iterative algorithm that alternates between identifying the phylo-genetic structure of the constituent clones comprising each tumour sample, anddetection of mutations leveraging the new phylogenetic structure.As sequencing costs continue to decrease (e.g. with Illumina’s NovoSeq plat-form), multi-sample whole genome sequencing of tumours will continue to prolif-erate (e.g. rapid autopsy program) as a viable experimental design. Thus, MuClonewill be an asset in the arsenal of analytical methods deployed to interpret evolu-tionary properties of cancer and to gain insights into clonal dynamics in time andspace.440.0 0.5 1.0False positive rate0.00.51.0True positive rateROC curve Sample P4a, AUC=1.0P4c, AUC=1.0P4e, AUC=1.0P4g, AUC=0.9875P4i, AUC=0.99913Figure 2.18: MuClone’s Roc curves and the area under the curve (AUC) forpatient 4. MuClone parameters are: wildtype prior = 0.5, ΦT = 0.02,tumour content = 0.75, error rate = 0.01, and precision parameter =1000.45MuCloneMultiSNVMutectNaive MuCloneStrelkaMutationSeqFreeBayes0.650.700.750.800.850.900.951.00PerformanceYouden’s IndexMuCloneMultiSNVMutectNaive MuCloneStrelkaMutationSeqFreeBayesSensitivityMuCloneMultiSNVMutectNaive MuCloneStrelkaMutationSeqFreeBayesSpecificity(a) (b) (c)Figure 2.19: Performance comparison of different methods on whole genomesequencing data for patient 7 with high grade serous ovarian cancer:(a) Youden’s index, (b) Sensitivity, and (c) Specificity across differ-ent mutation detection methods (from left to right: MuClone (darkblue), MultiSNV (orange), MuTect (light blue), Naive MuClone (yel-low), Strelka (purple), MutationSeq (brown), and FreeBayes (pink)).MuClone parameters are: wildtype prior = 0.5, ΦT = 0.02, tumourcontent = 0.75, error rate = 0.01, and precision parameter = 1000.460.0 0.5 1.0False positive rate0.00.51.0True positive rateROC curve Sample P3a, AUC=0.99917P3c, AUC=0.97942P3e, AUC=0.97799P3g, AUC=0.98231Figure 2.20: MuClone’s Roc curves and the area under the curve (AUC) forpatient 7. MuClone parameters are: wildtype prior = 0.5, ΦT = 0.02,tumour content = 0.75, error rate = 0.01, and precision parameter =1000.47MuCloneMultiSNVMutectNaive MuCloneStrelkaMutationSeqFreeBayes0.20.40.60.81.0PerformanceYouden’s IndexMuCloneMultiSNVMutectNaive MuCloneStrelkaMutationSeqFreeBayesSensitivityMuCloneMultiSNVMutectNaive MuCloneStrelkaMutationSeqFreeBayesSpecificity(a) (b) (c)Figure 2.21: Performance comparison of different methods on whole genomesequencing data for patient 9 with high grade serous ovarian cancer:(a) Youden’s index, (b) Sensitivity, and (c) Specificity across differ-ent mutation detection methods (from left to right: MuClone (darkblue), MultiSNV (orange), MuTect (light blue), Naive MuClone (yel-low), Strelka (purple), MutationSeq (brown), and FreeBayes (pink)).MuClone parameters are: wildtype prior = 0.5, ΦT = 0.02, tumourcontent = 0.75, error rate = 0.01, and precision parameter = 1000.48Figure 2.22: MuClone’s Roc curves and the area under the curve (AUC) forpatient 9. MuClone parameters are: wildtype prior = 0.5, ΦT = 0.02,tumour content = 0.75, error rate = 0.01, and precision parameter =1000.49MuCloneMultiSNVMutectNaive MuCloneStrelkaMutationSeqFreeBayes0.50.60.70.80.91.0PerformanceYouden’s IndexMuCloneMultiSNVMutectNaive MuCloneStrelkaMutationSeqFreeBayesSensitivityMuCloneMultiSNVMutectNaive MuCloneStrelkaMutationSeqFreeBayesSpecificity(a) (b) (c)Figure 2.23: Performance comparison of different methods on whole genomesequencing data for patient 10 with high grade serous ovarian cancer:(a) Youden’s index, (b) Sensitivity, and (c) Specificity across differ-ent mutation detection methods (from left to right: MuClone (darkblue), MultiSNV (orange), MuTect (light blue), Naive MuClone (yel-low), Strelka (purple), MutationSeq (brown), and FreeBayes (pink)).MuClone parameters are: wildtype prior = 0.5, ΦT = 0.02, tumourcontent = 0.75, error rate = 0.01, and precision parameter = 1000.50Figure 2.24: MuClone’s Roc curves and the area under the curve (AUC) forpatient 10. MuClone parameters are: wildtype prior = 0.5, ΦT = 0.02,tumour content = 0.75, error rate = 0.01, and precision parameter =1000.511 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19True cluster12345678910111213141516171819Predicted cluster0.000.040.080.120.16Figure 2.25: Classification of 153 mutations of patient 1 with high gradeserous ovarian cancer across 6 samples. Bin (i, j) shows the fractionof mutations in cluster i that were classified into cluster j by MuClone.MuClone parameters are: wildtype prior = 0.5, ΦT = 0.02, error rate= 0.01, tumour content = 0.75, and precision parameter = 1000. 93%of the elements are diagonal.52Total mutationsNumber of common mutations (a)(b)(c)Figure 2.26: Comparison of detected mutations from MultiSNV, MuTect,MuClone, and TRACERx on whole exome sequencing data fromnon-small cell lung cancer patients: CRUK0003 (dark blue),CRUK0004 (orange), CRUK0005 (dark purple), CRUK0013 (lightblue), CRUK0062 (light purple), CRUK0063 (brown), CRUK0065(pink), and CRUK0094 (grey). To illustrate overlap in detected mu-tations from all combinations of different methods, in panel (a) thenumber of mutations called in all selected methods, but not in anyother method, are plotted for each patient. Dark circles in the columnsof panel (c) indicate selected methods. For example, the first columnin panel (c) indicates mutations that are only called by MultiSNV (allcircles are grey except the one corresponding to MultiSNV). Panel (b)displays the total number of mutations called for each method. Mu-Clone removes reads with mapping quality less than 5 and positionswhich have (i) a variant nucleotide in a normal sample, (ii) more than40% filtered basecalls; or (iii) more than 75% of the reads that crossthe site have deletions in any of the samples. Other methods executedwith default settings. See UpSet tool for additional visualization de-tails [80].53Performance(a) (b) (c)Figure 2.27: MuClone’s performance with inaccurate clonal information:The mutations with full copy number information across all 8 non-small cell lung cancer patients are used as ground truth. To asses theperformance, we add noise from a zero mean normal with standard de-viation equal to 0 (dark blue), 0.01 (orange), 0.1 (light blue), and 0.2(purple) and also remove the clonal information of different fractionsof clusters. (a) Youden’s index, (b) Sensitivity, and (c) Specificity ofMuClone with different fraction of clusters information removed.54Chapter 3Single cell somatic mutationdetection through incorporationof the underlying phylogeny“Nothing truly valuable can be achieved except by the unselfishcooperation of many individuals.”- Albert EinsteinThis chapter introduces a new model, CellMutScope, to detect mutations (SNVs)in each cell based on the joint analysis of a large number (1000 - 10,000 or more)of whole genome single cell sequencing data. The model detects SNVs incorpo-rating the underlying CNV phylogenetic tree. Moreover, we expand the utility ofCorrupt model that infers the underlying phylogeny of a group of cells based onCNV data [1]. We frame how to incorporate SNV data in the Corrupt model (i) toinfer the underlying phylogeny using SNVs with or without CNV data, or (ii) to ex-tend the input underlying phylogenetic tree with the SNVs placed on the tree. Theextended tree discloses the cell-to-cell genomic variability information and showssimilar cells that are grouped together. Sections 3.1 and 3.2 define the problemand review the existing approaches, respectively. Section 3.3 reviews the fast scal-able Corrupt model that infers the underlying CNV tree. Section 3.4 frames how to55incorporate SNV data into the Corrupt model and introduces CellMutScope. Sec-tion 3.5 synthetically evaluates the results in terms of SNV calls and the phylogenywhen SNV data is also used.3.1 IntroductionSingle cell genomics is critical for deciphering the evolutionary process (clonal dy-namics) of cancer during cancer initiation, progression, and during ongoing evo-lution as a patient undergoes treatment. Bulk sequencing protocols sequence amixture of cells; computational methods are then applied to deconvolve the dataand to infer the underlying clones. However, it is challenging to resolve minorclones from bulk data. The challenge is because of the presence of sequencingerrors, and also because of sampling issues related to intra-tumour heterogeneityin solid tumours [141]. Single cell sequencing data promises to reveal the genomicdetails of single cells and the underlying tumour clones. Genomic variations (sin-gle nucleotide mutations, copy number variation, and structural variants) contain arecord of ancestral relationships between the cells, and the relationship is reflectedin a phylogenetic tree. Although single cell sequencing technologies provide adirect data to infer the phylogeny (to study the clones and their dynamics), ex-perimental challenges in capturing nuclei of individual cells and difficulties in se-quencing them with even coverage have delayed the large scale analysis of singlecell genomes [41, 75].The low quality of single cell sequencing data (SCS data) is rooted in the factthat a single cell only contains about 6-7 pg of genomic DNA while the typicalamount of DNA needed to construct a library in next generation sequencing (NGS)platforms is on the order of 1000 pg [21]. This means that the genome or partsof the genome needs to be amplified before sequencing [102, 109]. However,amplification of DNA is a stochastic process correlated with genome accessibil-ity and GC content [11]. In order to generate a sufficient quantity of DNA forsequencing, several protocols based on whole genome amplification (WGA) havebeen suggested [9, 40, 56, 75, 89, 128, 145]. While WGA permits sequencing tohigher coverage depth and breadth [37, 128, 141], it introduces coverage biases.In order to overcome the coverage bias issue, Zahn et al. have proposed a scalable56micro-fluidics based protocol, called Direct Library Preparation (DLP) [141]. Thisprotocol requires no pre-amplification and has lower biases compared to WGA pro-tocols. In the DLP protocol DNA of individual cells are first fragmented into shortsequencing inserts, and then through a few cycles of PCR, the sequencing adaptorsand index barcodes are added. This process produces exact duplicates that can beremoved computationally. The breadth and depth of DLP single cell sequencingdata coverage is lower than the bulk whole genome sequencing data. However, theresultant ‘bulk-equivalent’ data generated by aggregating cells has uniform cov-erage similar to the bulk data [141]. The uniform coverage accommodates CNVcalling, but the sparsity of the data introduces difficulties for SNV calling. In addi-tion, the isolation of cells may fail and multiple cells (mostly double cells) may besequenced together.The DLP data which is composed of thousands of cells, is sequenced shallowly(< 1x) across the genome. Following the uniform coverage of the DLP data, Zahnet al. proposed a tool for CNV calling [141]. Bouchard et al. then utilized the CNVdata to infer the underlying phylogeny in the Corrupt [1].In this chapter, we frame how to utilize SNV data in the Corrupt model (i) toinfer the underlying phylogeny, or (ii) to extend the underlying CNV tree by addingthe SNVs on the tree. The extended tree (that includes both CNVs and SNVs) showsthe genomic variability between individual cells and reveals cells that are groupedtogether. Incorporation of SNV data leads to a more distinguished identification ofminor clone structures, and this helps resolve the phylogeny in more detail. Theresultant phylogenetic tree contains a record of the ancestral relationships betweencells.The computational complexity of the Corrupt model tree inference at each iter-ation is linear with the number of cells and traits. Missing data in Bayesian modelsrequires a higher number of iterations before convergence [72]. This is also thecase for the above framework. We observe an undesirable increase in computa-tional time with missing rate in the SNV data. Therefore, it is reasonable to useonly the CNV data in tree inference step and add the SNVs afterwards.We propose a model, CellMutScope, to detect SNVs in every single cell. Weposit that incorporating the CNV data and the underlying phylogenetic tree acrossa large number of whole genome single cell sequencing data help overcome the57sparsity of single cell sequencing data. In contrast to MuClone’s output which is avector of length N (number of loci), CellMutScope output is a matrix of size M×N(M is the number of cells1. MuClone calls a mutation to be present if it is present inany of the cells in the sample, while CellMutScope calls the mutation status of eachcell individually. The overall consolidated framework that unites SNVs and CNVsover a phylogeny will improve the understanding of the tumour clonal dynamics.It also helps understanding of how the combination of different genetic variations(CNVs and SNVs) evolve to the observed results in the current population.3.2 Literature reviewMost of the existing methods utilize SCS data obtained from WGA protocols. Whilethese protocols are efficient in preparing enough material for the sequencing, theyhave high rate of allelic drop-out (about 10% -20%). The high allelic drop-out rateis the result of random amplification of only one allele at a heterozygous genotypesite. This complicates the SNV detection and results in a coverage bias. Moreover,existing methods are only focused on the analysis of SNVs across sites with enoughcoverage (roughly more than 5x).Monovar addressed the issue of low and uneven data coverage by poolingthe sequencing data across cells with an assumption that the sites are indepen-dent [140]. Monovar uses data with coverage between 6x to 24x, and it skips theones with less coverage. It has the asymptotic complexity of O(M3) for geno-typing M single cells. SCcaller identifies SNVs for each cell while accountingfor local allelic amplification biases [33]. The model requires the data coverageabout 5x, and it requires the single nucleotide polymorphisms (SNPs) data. Thisdata is not always available. Furthermore, it can recover mutations neither fromdrop-out events nor from loss of heterozygosity [33]. Conbase and LiRA utilizeread-phasing to correct for errors and allelic drop-outs [12, 54]. Conbase unlikeLiRA performs joint variant calling across the population of cells. The use of readphasing decreases false discovery rate and increases specificity. However, suchanalysis is only possible for regions in proximity to SNPs. This prevents callingSNVs across all of the bases covered by the reads. SCIΦ simultaneously detects1M represents the number of samples in Chapter 2. A sample in this chapter represents a cell.58(a) (b) Figure 3.1: (a) A schematic view of CNV tree where the black nodes representcells, and the blank nodes represent CNV traits. (b) A schematic view ofthe tree after adding SNV loci. The grey nodes represent SNV traits.SNVs and infers the underlying phylogeny. Although SCIΦ claims to detect muta-tions in single cells with very low or even no coverage, it requires at least two cellsto show an alternative nucleotide count of at least three. The runtime complexity ofSCIΦ is O(X ×max(MN,Cu)), where Cu is the number of unique coverage valuesof the experiment for M cells and N loci, and X is the number of iterations [115].We propose a method, CellMutScope, that calls SNVs given the underlyingphylogeny. Its computational complexity is O(M+N). With this scalable frame-work each SNV locus is quickly processed. Inferring the tree from CNV data itselfhas a computational complexity of O(M +N) per iteration [1]. Therefore, theunified framework provides a comprehensive genomic data analysis tool for largescale low-coverage genome-wide single cell data.Figure 3.1(a) shows a schematic view of a CNV tree in which the black nodesrepresent cells, and white nodes represent CNV traits. Figure 3.1(b) shows the treeafter running CellMutScope and adding SNV loci (The grey nodes represent SNVtraits).593.3 Corrupt modelMost of the copy number analysis methods start with dividing the genome intonon-overlapping segments. Each segment is assumed to be homogeneous in copynumber, and the number of overlapping reads is proportional to the DNA copynumber at the corresponding segment. The issue of coverage bias in SCS data ob-tained from WGA, hinders development of methods to utilize CNV data. Navinet al. and Wang et al. infer the copy number of single cells from WGA based SCSdata [89, 128]. The DLP data is more uniform and has lower coverage bias com-pared to WGA protocols. Therefore, the data is more suitable for detecting copynumber alterations. Zahn et al. developed a method that exploits DLP data to inferthe copy numbers [141]. Following this work, Bouchard et al. have proposed Cor-rupt model that incorporates CNV data and infers the underlying phylogeny withO(M+N) per iteration [1].Corrupt is a Bayesian framework for inferring phylogenetic trees using copynumber information from large scale low-coverage genome-wide single cell data.A schematic view of the tree is depicted in Figure 3.1(a). In order to simplify thesite dependencies, Corrupt encodes CNV data into binary change points. The copynumber data is indexed by genomic bins, and each bin has an integer copy numbervalue. In the binary change point representation of data, the data is indexed bythe ‘space’ between two adjacent bins; and it is a change point if there is a copynumber change between two bins. Let us define n as a locus between two adjacentbins, and m as an index for the cells. The copy number value of the bins beforeand after locus n in cell m is denoted by c−m,n and c+m,n, respectively. Let M and Ndenote the disjoint set of cells and loci, respectively. The trait value for cell m ∈Mand locus n ∈ N is defined as:xm,n = 1[c−m,n 6= c+m,n] ={1, if c−m,n 6= c+m,n,0, otherwise.(3.1)The matrix x = (xm,n) is not observed directly, since c−m,n, c+m,n are unobserved.Therefore, the proposed observation probability model is p(y|x,θCNV ), where yis the observed data and θCNV encodes the error rates in copy number calls. Letus define θCNVn = (rFPm,n(θ),rFNm,n(θ)) to represent the locus specific false positive60and false negative rates of cell m. The error rates are rFPm,n(θ) ∈ (0,1) and rFNm,n(θ) ∈(0,1). For simplicity, the likelihood probability distribution denoted by p(yCNVm,n |xCNVm,n ,θCNV )is written as:p(ym,n|xm,n,θ) = p(ym,n|xm,n,rFPm,n(θ),rFNm,n(θ)) = erFPn (θ),rFNn (θ)xm,n,ym,n , (3.2)where the error matrix iserFPm,n(θ),rFNm,n(θ) =[1− rFPm,n(θ) rFPm,n(θ)rFNm,n(θ) 1− rFNm,n(θ)].In the simplest case, these parameters are globally constant across all cells and loci,whererFPm,n(θ) = rFPrFNm,n(θ) = rFN .Therefore,p(y|x,θ) = ∏m∈M∏n∈Np(ym,n|xm,n,rFP(θ),rFN(θ)), (3.3)= (rFP)nFP(rFN)nFN(1− rFP)nN−nFN (1− rFN)nP−nFP , (3.4)where nFP and nFN denotes the number of false positive and false negative in-stances, respectively. They are computed as:nFP = ∑m∈M∑n∈NI[xm,n = 0,ym,n = 1],nFN = ∑m∈M∑n∈NI[xm,n = 1,ym,n = 0],where, nP and nN are the number of positive and negative instances, withnp = ∑m∈M∑n∈NI[ym,n = 1],61andnN = |M||N|−nP.Initially, the parameters are estimated using the whole matrix, then at each iteration,the parameters are updated by comparing the row or the column that is modifiedcompared to the previous iteration (each iteration or move modifies only one row orcolumn of the matrix x). In the Corrupt inference model, the error rate parameterscan be set to local or global.3.3.1 Probability modelAn underlying assumption in the Corrupt is the validity of perfect phylogeny in thelatent data (xm,n), and not in the observed data (ym,n). The model can be describedby a two-step generative process: (i) sampling a mutation tree, and (ii) samplingcell assignments.Sampling a mutation tree Let T denote a set of tz trees spanning node set V ,where V includes root node v∗ plus a node for each of the N loci. The root nodev∗ puts an implicit direction on the edges of the tree. There is a directed path fromnode/locus n to n′ in tz, if and only if the trait indexed by n is emerged in a cellprior to the trait indexed by n′.Sampling cell assignments Assign each cell to a node in tz. If cell m is as-signed to locus n, then the cell has all the traits on the shortest path from locus n toroot v∗. If a cell is assigned to node v∗, then it doesn’t have any of the traits on thetree.3.3.2 InferenceThe goal is to explore the space of trees. Corrupt model proposes a move withcomputational complexity of O(N+M) to explore the exponentially large numberof neighbours. The move consists of four steps:1. Remove a locus (trait) node n ∈ N from the tree.2. Select a node v that is not a cell node on the tree from which node n isremoved.3. Pick a (possibly empty) subset from the children of v.62w v a b c d e a b c d e v Edge insertion Edge removal Figure 3.2: A schematic of removing(adding) a node from(to) the tree. Thebold edge and node w are removed, and the nodes connected to w isconnected to the parent node v. The triangles represent a subtree rootedat an specific internal node. [modified figure from [1]]4. Add a node w to represent locus n by adding an edge from v to w. Move thesubset of selected nodes from previous step to hang from w (See Figure 3.2.).The above steps are followed when sampling or maximizing. In order to select anode, the following probability is efficiently calculated:ρ¯v =ρv∑v˜∈Rρv˜, (3.5)where v˜ ∈ R, and R = {v∗}⋃N\{n}. The value of ρv isρv = ∑t∈N nv (t\n)p(t)p(y|x(t),θ), (3.6)where N nv(t\n)denotes the set of neighbour trees of the tree t with locus node nremoved. After further simplification, ρ¯v is estimated asρ¯v =ρv∑v˜∈Rρv˜,=(∏vi∈Children(v)(p0vi+p1vi)p0v)∑v¯∈R(∏v¯i∈Children(v¯)(p0v¯i+p1v¯i)p0v¯) , (3.7)63where for all other v ∈ R and s ∈ {0,1},psv = ∏v′′∈Children(v)psv′′ . (3.8)For all nodes m corresponding to a cell and s ∈ {0,1}, the likelihood probabilityfunction is defined as:psm,n = p(ym,n|xm,n,s,rFPn (θ),rFNn (θ)). (3.9)After the attachment node v is selected, for each child vi of v, a Bernouli dis-tribution with success parameter p1vip1vi+p0viis considered. The success corresponds tomoving vi into a child of the newly re-introduced node w. The realization of theBernoulli variables can be sampled from for approximating the posterior distribu-tion. Also, it can be maximized for adding a node to the tree.3.3.3 Approximation of the posterior distributionThe posterior distribution ispi(t,θ) ∝ p(t)p(θ)p(y|x(t),θ).Here the deterministic function is denoted by x(t). Given tz ∈ T , the matrix X isa deterministic function obtained by setting xm,n = 1 if node m is a descendant ofnode n in tz, and zero otherwise. Different methods are used to approximate theposterior function:1. A direct Markov chain Monte Carlo (MCMC) scheme.2. A non-reversible Parallel Tempering (PT) algorithm [133].3. The sequential change of measure (SCM) scheme from [28] with the adap-tive scheme from [144].643.3.4 Summary via minimum Bayes estimatorMinimum Bayes risk estimator is used to summarize the posterior distribution:argmint∈T∑t ′∈T∫L(t, t ′)pi(t, dθ), (3.10)where t ′ ∈T , andL(t, t ′) = ∑n∈N∑m∈M|xm,n(t)− xm,n(t ′)|.This minimization problem itself can be solved through a greedy procedure. Moredetails can be found in [1].3.4 MethodThe Corrupt model utilizes CNV data for inferring the underlying phylogeny. InSection 3.4.1 we frame how to use SNV data (allelic read counts) with or withoutCNV data (i) to infer the underlying phylogeny, or (ii) to add SNVs on the un-derlying CNV tree. With this framework, we can infer the underlying phylogenytree using both CNV and SNV data. This latter approach (adding the SNV loci to thetree) helps skip the computationally expensive inference step (when both CNVs andSNVs are used in tree inference) and infers the extended tree (with both CNVs andSNVs) in a time efficient manner. Next, we propose a new model, CellMutScope,for detecting SNVs in every single cell through adding the SNV loci to the under-lying tree. Incorporating the underlying phylogeny and CNV data help overcomethe high missing rate of the data and detects SNVs in every individual cell. Theresultant extended tree and the SNV calls discloses further detail about the clonalstructure of the cells that were not distinguishable based on only CNV data.3.4.1 Incorporating single nucleotide data in the Corrupt modelThe Corrupt model takes binary representation of copy number data. This datadepends on the unobserved copy number variables to infer the phylogeny. Theproposed observation probability model is defined in Equation 3.8. In this section,we posit an observation probability model for SNV data to extend the utility of the65Corrupt model, where the observed data is ySNVm,n = (bm,n,dm,n,cm,n). For locus nin cell m, the values of dm,n, bm,n and cm,n respectively represent: the total numberof reads covering the locus, the number of reads with a variant allele (comparedto the reference genome), and copy number. The observed data, in particular bm,ndepends on the mutation status of cell m at locus n, denoted by xSNVm,n . The errormodel is θ SNV = (εFP,εFN), where εFP and εFN are false positive rate and falsenegative rate, respectively. Thereforeqsm,n = p(ySNVm,n |xSNVm,n ,θ SNV ) = p(bm,n|dm,n,cm,n,xSNVm,n = s,θ SNV ), (3.11)where dm,n and cm,n are given. The likelihood probability of cell node m is denotedby qsm,n, where s ∈ {0,1}. For s = 1, qsm,n reflects the likelihood of cell m beingmutated at locus n; and for s = 0, qsm,n reflects the likelihood of cell m not beingmutated at locus n.The probability qsm,n follows a mixture of binomial distributions depending onall possible genotype states of locus n at cell m. Given the copy number cm,n,the possible genotype states are G = {A . . .A,AA . . .B,A . . .BB, . . . ,B . . .B}, whereeach element has a length equal to cm,n. For example, the genotype AAB refers to agenotype with one variant allele B and 2 reference allele A. For each genotype stategi, where i indexes the elements of G , the mean parameter of the correspondingbinomial distribution is denoted by ξ im,n:ξ im,n =B(gi)cm,n, 1≤B(gi)< cm,n,1− εFP, B(gi) = cm,n,εFP, otherwise,(3.12)whereB(gi) represents the number of variant allele of genotype gi. Therefore, fors = 1,q1m,n = p(bm,n|dm,n,cm,n,xSNVm,n = 1,θ SNV ) (3.13)=cm,n∑i=1p(gi)[ξbm,nm,n (1−ξm,n)dm,n−bm,n ] (3.14)+ εFN [εbm,nFP (1− εFP)dm,n−bm,n ].66The value of p(gi) equals 1−εFNcm,n , and εFN represents the error due to mutation lossor tree errors.If the mutation status of cell m at locus n is wildtype (i.e mutation is notpresent), then the possible genotype states should not have any variant allele. Theonly possible genotype state is {A . . .A}. The mean parameter of the binomialdistribution equals εFP (false positive rate). Therefore,q0m,n = p(bm,n|dm,n,cm,n,xSNVm,n = 0,εFP). (3.15)Algorithm 1 summarizes the assumed model to generate the data where gm,n is thetrue genotype of locus n of cell m.With the proposed probability model for SNVs, we can incorporate both SNVdata and CNV data to infer the underlying tree phylogeny in the Corrupt model.Therefore, Equation 3.4 is updated as:p(y|x,θ) = ∏m∈M∏n∈NCNVp(yCNVm,n |xCNVm,n ,θCNV ) ∏n∈NSNVp(ySNVm,n |xSNVm,n ,θ SNV ), (3.16)where M and N is the disjoint set of cells and loci, respectively. N includes bothCNV and SNV traits, and N = NCNV +NSNV . Therefore, Equation 3.7 (when locusn is removed from the tree) can be rewritten in the following form:ρ¯v =(∏vi∈Children(v)(γ0vi+γ1vi)γ0v)∑v¯∈R(∏v¯i∈Children(v¯)(γ0v¯i+γ1v¯i)γ0v¯) , (3.17)where γsv , for s ∈ {0,1} is:γsv ={psv, If n represents a CNV loci,qsv, If n represents a SNV loci.67Algorithm 1: Data’s underlying generative modelInput: dm,n,cm,n,θ SNV ,s;Output: bm,n;if s = 1 thendraw l ∼ Bernoulli(εFN);if l = 0 thendrawB(gm,n)∼ Uniform([1,cm,n]) ;ifB(gm,n) = cm,n thendraw bm,n ∼ Binomial(dm,n,1− εFP);elseξm,n =gm,ncm,n;draw bm,n ∼ Binomial(dm,n,ξm,n);endelsedraw bm,n ∼ Binomial(dm,n,εFP);endelsedraw bm,n ∼ Binomial(dm,n,εFP);endreturn bm,n;For v ∈ R = {v∗}⋃N\{n}, and s ∈ {0,1} , the value of qsv isqsv = ∏v′′∈Children(v)qsv′′ . (3.18)For the cell nodes that are the leaves of the tree qsv = qsm,n.3.4.2 Detection of SNVs at every individual cellGiven the underlying tree (denoted by t) and the read counts data (denoted by y),here the goal is to calculate the posterior probability of xm,n, where xm,n denotesthe mutation status of locus n at cell m. According to the underlying tree, cells areattached to the trait nodes, and they are leaves of the tree (See Figure 3.1 (b)). The68joint probability distribution of xm,n, y and t can be written as:p(xm,n,y, t) = ∑v∈R∑t ′∈N nv (t\n)p(xm,n, t ′,y) (3.19)= ∑v∈R∑t ′∈N nv (t\n)p(xm,n|t ′)p(y|t ′)p(t ′), (3.20)where R is the set of all loci nodes in the tree (including the root) excluding locusn. The joint probability distribution is calculated asp(xm,n = 1,y, t) = ∑v∈P(m,t)∑t ′∈N nv (t\n)p(y|t ′)p(t ′), (3.21)where xm,n is a deterministic function of t ′ which belongs to N nv (t\n). The valueof xm,n is 1 if locus n is an ancestor of cell m, otherwise it is 0. The set P(m, t)denotes all the nodes on the path from cell m to the root of the tree (including theroot and excluding the cell m node). An example of the path on an imaginary tree isdepicted in Figure 3.3. The nodes coloured in green belong toP(m, t). Therefore,the posterior probability distribution of xm,n = 1 yieldsp(xm,n = 1|y, t) = p(xm,n = 1,y, t)p(y, t) =∑v∈P(m,t)∑t ′∈N nv (t\n) p(y|t ′)p(t ′)p(y, t). (3.22)Rewriting Equation 3.22 assuming uniform probability distribution for p(t ′) yields:p(xm,n = 1|y, t) ∝ ∑v∈P(m,t)∑t ′∈N nv (t\n)p(y|t ′),= ∑v∈P(m,t)∑t ′∈N nv (t\n)∏n′∈N∏m′∈Mp(ym′,n′ |t ′),= ∑v∈P(m,t)∑t ′∈N nv (t\n)∏n′∈Nn′ 6=n∏m′∈Mp(ym′,n′ |t ′) ∏m′∈Mp(ym′,n|t ′),= K1 ∑v∈P(m,t)∑t ′∈N nv (t\n)∏m′∈Mp(ym′,n|t ′),= K1 ∑v∈P(m,t)∑t ′∈N nv (t\n)∏m′∈M\vp(ym′,n|t ′) ∏m′∈Mvp(ym′,n|t ′),69where N denotes the set of all trait nodes, M denotes the set of all cell nodes, Mvdenotes the cells that are a descendant of node v, and M\v denotes the cells thatare a not descendant of node v. The product of the likelihood contributions fornon-descendant nodes can be calculated by taking the product of q0m for all cells,divided by the ones that are descendant of v:∏m′∈M\vq0m′ =q0v∗q0v.Therefore:p(xm,n = 1|y, t) ∝ K1 ∑v∈P(m,t)q0v∗q0v∑t ′∈N nv (t\n)∏m′∈Mvp(ym′,n|t ′). (3.23)The likelihood contribution of descendant cells can be re-indexed by a binaryvector of s = (s1,s2, . . . ,sk), where si ∈ {0,1}, and si = 1 if the child v is to bemoved into a child of the node n. The value of k denotes the number of children ofv. The i∗th child of v which is on the path from node v to cell m is called v∗i . Thisimplies si∗ = 1 (See Figure 3.3). Therefore:∑t ′∈N nv (t\n)∏m′∈Mvp(ym′,n|t ′) = q1v∗m1∑s1=01∑s2=0. . .1∑si−1=01∑si+1=0. . .1∑sk=0k∏i=1i6=i∗qsivi . (3.24)Rewriting Equation 3.23 using Equation 3.24 yields:p(xm,n = 1|y, t) ∝ K1 ∑v∈P(m,t)q0v∗q0vq1v∗m1∑s1=01∑s2=0. . .1∑si−1=01∑si+1=0. . .1∑sk=0k∏i=1i 6=i∗qsivi ,= K1 ∑v∈P(m,t)q0v∗q0vq1v∗mk∏i=1i 6=i∗(q0vi +q1vi),= K1 ∑v∈P(m,t)q0v∗q0v∏ki=1(q0vi +q1vi)(q0vi∗ +q1vi∗ )q1vi∗ ,= K1q0v∗ ∑v∈P(m,t)q1vi∗q0v(q0vi∗ +q1vi∗ )k∏i=1(q0vi +q1vi). (3.25)70m Root v v1 v2 vi* vk Figure 3.3: A schematic view of the underlying tree inferred from CNV andSNV loci across multiple cells. Black and white nodes represent cellsand loci, respectively. The grey triangle represents a subtree rooted at anode. It includes all of the nodes and edges in the subtree.71With this approach, the posterior probability of xm,n is calculated based onplacing locus n at different places on the tree. This considers different subset ofcells (sub-trees of the underlying tree) to calculate the posterior probability of xm,n.3.4.3 Model’s asymptotic complexityThe computational complexity of 3.25 is O(MN) with M number of cells and Nnumber of loci. In order to reduce the complexity of calculating p(xm,n = 1|y, t) foreach locus per cell,P ′(m, t) is defined to denote the nodes sitting on the path fromroot to cell m, excluding the root node and including the cell m node. Then,q∗v =k∏i=1(q0vi +q1vi). (3.26)Therefore,K1q0v∗ ∑v∈P(m,t)q1vi∗q0v(q0vi∗ +q1vi∗ )k∏i=1(q0vi +q1vi) = K1q0v∗ ∑v∈P ′(m,t)q1v(q0v +q1v)q∗parent(v)q0parent(v).Calculating p(xm,n = 1|y, t) with a recursive approach reduces the complexityfrom O(MN) to O(M+N), where N = nSNV +nCNV .3.5 Benchmarking experiments3.5.1 Synthetic data generationTo generate the synthetic data, a tree is uniformly sampled from the set of undi-rected trees with N = NSNV +NCNV number of nodes. Then cells are assigned tothe trait nodes such that the cell nodes are the leaves of the the tree. The internalnodes are randomly assigned to either a SNV locus or to a CNV locus. For eachSNV locus, the copy number and the number of reads covering the locus (depth)are independently generated. The copy number is generated from a discrete uni-form distribution from 0 to cmax (cmax equals 5 in the analysis); the depth isgenerated from a Poisson distribution given the coverage. Given the depth and thecopy number, a random genotype is uniformly sampled from all possible genotype72states (G ). The number of variant reads is sampled from a binomial distributiongiven depth, copy number, false positive and false negative error rate parameters(See Algorithm 1).3.5.2 Synthetic data evaluationIn this section we evaluate both the topology of the inferred tree, and the accuracyof SNV calls. We define a distance D to evaluate the tree topology performance.Consider X to be a matrix derived from the phylogeny, where x(i, j) = 1 if locus iis an ancestor of cell j. If X1 is derived based on a reference tree and X2 is derivedfrom an inferred tree, then the distance between X1 and X2 is estimated byD =∑i∑j|(x1(i, j)− x2(i, j))|/z, (3.27)where z is the total number of entries of matrix X1 or X2 (both are the same size).The performance of SNV calls are estimated by sensitivity, specificity andYouden’s index. They are estimated given both the binary SNV calls, and theground truth calls (derived from the reference tree). The SNV calls are binarized bysetting a threshold equal to 0.5 (probability greater than 0.5 is 1 and smaller than0.5 is 0). The topology of the reference tree provides the ground truth of SNV callsin which locus i is mutated in the entire descendant cells of a sub-tree rooted at i.In Figure 3.4 and 3.5, we demonstrate the performance of SNV calls for datawith different coverages. The data is generated for 200 SNV loci, 100 CNV lociacross 1000 cells. We randomly picked 90% of SNV loci and utilized them withthe CNV loci for inferring the tree. The remaining 10% SNV loci are added after-wards. The nodes are assigned to 10 random tree topologies. We assign a copynumber to each SNV locus. The copy number is uniformly sampled from {1 . . .5}.The genotype is also uniformly sampled from all possible genotype states giventhe copy number. The number of reads and the number of variant alleles are gener-ated following the binomial probability distribution defined in 3.15 and 3.13 giventhe genotype and different coverage means of 0.04, 0.07, 0.1, 0.5 or 2. The in-ference is done for 10 different seeds, and non-reversible Parallel Tempering (PT)algorithm [133] with 1000 and 5000 number of iterations with one chain [1]. The73●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ● ● ●Sensitivity Specificity0.04 0.07 0.1 0.5 2 0.04 0.07 0.1 0.5 20.000.250.500.751.000.000.250.500.751.00CoverageMethodMCMC−PTMCMC−PT−5kCoverageFigure 3.4: The performance (Sensitivity and Specificity) of SNV calls acrossdata coverages of 0.04,0.07,0.1,0.5, and 2. The data is generated for200 SNV loci, 100 CNV loci across 1000 cells for 10 different topolo-gies. Only 90% of SNV loci are used for the inference, and the remainingloci are added to the tree afterwards. Copy number is uniformly sam-pled from {1 . . .5}. MCMC-PT and MCMC-PT-5k denote the inferencealgorithm using PT with 1000 and 5000 iterations, respectively.values of εFP and εFN are equal to 0.01. The inference parameters are constantthrough all the following analysis. As expected, the higher coverage the betterperformance. As depicted in Figure 3.4 and 3.5, the performance is not very highwith 0.07x coverage. However, the performance reaches almost 1 with 2x cover-age. Note that the results are evaluated per cell, and these coverages are all lowerthan the coverage of the data used in all other existing methods. Figure 3.6 de-picts the performance of the model (in terms of tree topology) improves with morecoverage.Next, we investigate the performance of the CellMutScope when different por-tion of SNV loci are used for inferring the tree and the remaining loci are added tothe inferred tree. Using a subset of SNVs is helpful for reducing the run time. Thedata is generated for 200 SNV loci, 100 CNV loci across 1000 cells. The data cov-erage is 0.07 and cmax is equal to 5. The percentage of loci used for the inferenceare 20%, 50% and 90%. In Figure 3.7 and 3.8, we explore how the performancechanges with different proportion of SNVs used for tree inference. The perfor-mance of detecting SNV is higher when 90% percent of SNV loci are provided.However, the difference is not significant. Figure 3.9, suggests that incorporatingSNVs improves resolving the phylogeny. However, for larger number of the nodes,we may need to increase the number of iterations to reach a higher accuracy.74●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●● ●Youden’s index0.04 0.07 0.1 0.5 2−0.50.00.51.0CoverageMethodMCMC−PTMCMC−PT−5kFigure 3.5: The performance (Youden’s index) of SNV calls across data cov-erages of 0.04,0.07,0.1,0.5, and 2. The data is generated for 200 SNVloci, 100 CNV loci across 1000 cells for 10 different topology. Only 90%of SNV loci are used for the inference, and the remaining are added tothe tree afterwards. Copy number is uniformly sampled from {1 . . .5}.MCMC-PT and MCMC-PT-5k denote the inference algorithm using PTwith 1000 and 5000 iterations, respectively.●●●●●●●●●●●● ●●D0.04 0.07 0.1 0.5 20.010.020.030.040.05CoverageMethodMCMC−PTMCMC−PT−5kFigure 3.6: The performance (D) of tree topology across data coverages of0.04,0.07,0.1,0.5, and 2. The data is generated for 200 SNV loci, 100CNV loci across 1000 cells for 10 different topology. Only 90% of SNVloci are used for the inference, and the remaining are added to the treeafterwards. Copy number is uniformly sampled from {1 . . .5}. MCMC-PT and MCMC-PT-5k denote the inference algorithm using PT with1000 and 5000 iterations, respectively.75●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●Sensitivity Specificity0.2 0.5 0.9 0.2 0.5 0.90.250.500.751.000.20.40.6Precentage of SNV LociMethodMCMC−PTMCMC−PT−5kPrecentage of SNV lociFigure 3.7: The performance (Sensitivity and Specificity) of SNV calls if vari-ous percentage of SNV loci are used for the inference, and the remainingare added afterwards. The data is generated for 200 SNV loci, 100 CNVloci across 1000 cells for 10 different topology. Copy number is uni-formly sampled from {1 . . .5} and the coverage mean is 0.07. MCMC-PT and MCMC-PT-5k denote the inference algorithms, and use PT with1000 and 5000 iterations, respectively.●●●●●● ●●●● ●●●●●●●●●●●●●●●●●●●●Youden’s index0.2 0.5 0.9−0.50−0.250.000.250.500.75MethodMCMC−PTMCMC−PT−5kPercentage of SNV lociFigure 3.8: The performance (Youden’s index) of SNV calls if various per-centage of SNV loci are used for the inference, and the remaining areadded afterwards. The data is generated for 200 SNV loci, 100 CNV lociacross 1000 cells for 10 different topology. Copy number is uniformlysampled from {1 . . .5} and the coverage mean is 0.07. MCMC-PT andMCMC-PT-5k denote the inference algorithm, and use PT with 1000and 5000 iterations, respectively.76D0.2 0.5 0.90.0300.0350.0400.0450.050MethodMCMC−PTMCMC−PT−5kPercentage of SNV lociFigure 3.9: The performance (D) of tree topology if various percentage ofSNV loci are used for the inference, and the remaining are added af-terwards. The data is generated for 200 SNV loci, 100 CNV loci across1000 cells for 10 different topology. Copy number is uniformly sampledfrom {1 . . .5} and the coverage mean is 0.07. MCMC-PT and MCMC-PT-5k denote the inference algorithms, and use PT with 1000 and 5000iterations, respectively.In Figures 3.10 and 3.11, we explore how the number of cells affects the per-formance of SNV calls. The data is generated for 200 SNV loci, 100 CNV lociacross 100, 200, 500, 1000, 5000 and 10000 cells. The coverage mean is 0.07 andonly 20% of SNV loci are used for the tree inference. The other parameters staysimilar to the previous analysis. For this number of loci, there is a slight increasein the performance by increasing the number of cells from 100 to 1000. However,further increase of the number of cells slightly drops the performance. This canbe explained by the need to increase the number of iterations for tree inferencewith larger number of cells. Since the SNV data has a lot of missing values, thisnegatively affects the inference [72].3.6 ConclusionWe studied the use of tree phylogeny for the purpose of SNV detection across alarge number of clonally related cells. The proposed statistical framework uses thecopy number profile of each cell and the underlying phylogeny inferred based onthe copy number data. Our proposed method, CellMutScope, is fast and scalable77●●●●●●●●●●●●● ●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●Sensitivity Specificity100 200 500 1000 5000 10000 100 200 500 1000 5000 100000.250.500.751.000.00.20.40.6Number of CellsMethodMCMC−PTMCMC−PT−5kNumber of CellsFigure 3.10: The performance (Sensitivity and Specificity) of SNV callsacross different number of cells. The data is generated for 200 SNVloci, 100 CNV loci across 100, 200, 500, 1000, 5000, and 10000 cellsfor 10 different topology. Copy number is uniformly sampled from{1 . . .5}, and the coverage mean is 0.07. Only 20% of SNV loci areused for tree inference. MCMC-PT and MCMC-PT-5k denote the in-ference algorithms, and use PT with 1000 and 5000 iterations respec-tively.to the large number of cells. With the benefit of placing SNVs (when SNVs areincorporated into the Corrupt model) on the underlying tree beside CNV traits, itenables us to detect minor clones that were not distinguishable with the CNV data.The model uses DLP data which has lower depth coverage compared to the dataused in other existing methods. In order to increase the coverage for the data,computationally similar cells can be grouped together such that the pseduo-cellshave higher coverage (cells grouped together called as a pseduo-cell). Anotherapproach would be an experimental increase in the number of lanes at the cost of adecrease in the number of cells (that keeps experimental cost fixed.).Detection of SNVs at single cell level facilitates identification of clone specificdriver mutations using transcriptomic single cell sequencing data. As the field ma-tures, we expect that the method presented here to be incorporated into analyticallycomprehensive modelling of single cell whole genome sequencing data (speciallywhen thousands of cells are used to infer the tumour dynamics). For future work wesuggest incorporating other genomic aberrations like indels and structural variantsin reconstructing the underlying tree phylogeny.As sequencing cost continues to decrease and the quality of single cell sequenc-ing data continues to increase, single cell sequencing of multiple clonally related78●●●●●●●●● ●●●●●●●●●●●●●Youden’s index100 200 500 1000 5000 10000−0.50−0.250.000.250.500.75MethodMCMC−PTMCMC−PT−5kNumber of CellsFigure 3.11: The performance (Youden’s index) of SNV calls across differentnumber of cells. The data is generated for 200 SNV loci, 100 CNVloci across 100, 200, 500, 1000, 5000, and 10000 cells for 10 differenttopology. Copy number is uniformly sampled from {1 . . .5}, and thecoverage mean is 0.07. Only 20% of SNV loci are used for tree infer-ence. MCMC-PT and MCMC-PT-5k denote the inference algorithms,and use PT with 1000 and 5000 iterations, respectively.cells will proliferate as a valuable experimental design. Thus CellMutScope willbe an asset in the toolbox of analytical methods exploited to interpret tumour dy-namics in time and space and to enhance the insight into tumours evolutionaryprocess.79●●●D100 200 500 1000 5000 100000.030.040.050.060.07MethodMCMC−PTMCMC−PT−5kNumber of CellsFigure 3.12: The performance (D) of tree topology across different number ofcells. The data is generated for 200 SNV loci, 100 CNV loci across 100,200, 500, 1000, 5000, and 10000 cells for 10 different topology. Copynumber is uniformly sampled from {1 . . .5}, and the coverage mean is0.07. Only 20% of SNV loci are used for tree inference. MCMC-PTand MCMC-PT-5k denote the inference algorithms, and use PT with1000 and 5000 iterations, respectively.80Chapter 4Real data application“Each piece, or part, of the whole of nature is always merely anapproximation to the complete truth, or the complete truth so far aswe know it. In fact, everything we know is only some kind of approxi-mation because we know that we do not know all the laws as yet.”- Richard P. FeynmanThis chapter presents the results of analysis base on SNV data using singlecancerous cell sequencing data. Section 4.1 discusses the impact of accurate SNVdetection on deciphering the clonal dynamics of tumours. Section 4.2 briefly de-scribes the real data analysed in this chapter. Section 4.3 explains the initial datacleaning and main analysis. It also suggests a quality control step to improve thequality of the results. Finally, Section 4.4 shows the results and the genomic inter-pretations.4.1 IntroductionCancer is driven through the accumulation of genetic mutations at cellular level.Mutations can disrupt the regular activities of proteins. They can lead to a higherproliferation rate, if they confer selective advantage. A specific mutation may dif-ferently affect various tumours and patients. The difference in their impact leads tointra- and inter- heterogeneity of tumour cells, and that can explain the differencein therapy responses and survival rates.81Sequencing technologies provide data toward the goal of revealing the associa-tions between mutations and clinical outcomes. Insight on tumour clonal dynamicsoffers a new perspective to improve the diagnosis, prognosis and treatment deci-sions beyond conventional approaches [93].Mutations that occur in high impact genes or non-synonymous coding genescan be the underlying cause of cells phenotypes. Several genes such as TP53,BRCA1, BRCA2, PTEN, pRb, etc. are well known tumour suppressor genes whichinhibit cell proliferation and tumour development [5, 97]. Another group of genes,reported as housekeeping genes, like ACTB, GAPDH, AHSP, B2M, etc. governor prevent cell growth [97, 114]. Mutations that affect the regular activity of thesegenes promote unconditional cell growth. For example, recent studies reveal thedriving role of genes PIK3CA, PTEN, AKT1, PIK3CA, ERBB2, FGFR1, MYC,MAP3K1, TP53, CCND1 and RB1 in breast cancer [5, 65, 132], CDKN2A, KRAS,PTEN, RB1, CCNE1, EVI2A and LCP2 in ovarian cancer [45] and SPOP, FOXA1and MED12 in prostate cancer [7].In this chapter, we focus on the analysis of three real data sets, each com-promising thousands of clonally related cancerous cells. The model imputes themissing values of the data and reveals that clone specific mutations occurred inhigh impact genes. We evaluate the performance and report the high impact genesidentified in these data sets.4.2 DataThe first sample is from a patient with a high-grade serous ovarian cancer. Thesamples are across different points in time and in anatomical spaces. The othertwo data sets are from patients with breast cancer tissues xenografted in mice andtransplanted across multiple time points.4.2.1 Ovarian cancerHigh-grade serous ovarian cancer that arises from the serous epithelial layer inthe abdominopelvic cavity is the most aggressive subtype of ovarian cancer [13].Tissue samples from female patients with high-grade serous ovarian cancer areused to sequence cell lines in more than 1966 libraries [76]. The cell lines are from82the 3 clonally related high grade serous samples, sourced from one primary tumour(SA1090) and two relapse specimens after therapy (SA921 and SA922). SA1090 (orcalled OV 2295) and SA922 (or called OV 2295(R)) are from ascites, and SA921 (orcalled TOV 2295) is from the solid tumour. For simplicity, these three samplestogether are called OV 2295.4.2.2 Breast cancer xenograft samplesBreast cancer is the most common cause of death in women diagnosed with can-cer [110]. Medical practitioners increasingly use hormone receptors to categorizebreast cancer. Hormone receptors of healthy breast cells receive messages and pro-vide instructions to cells for their growth and functionality [2]. A cancer is calledestrogen-receptor-positive (or ER+), if it has receptors for estrogen. It is calledprogesterone-receptor-positive (PR+), if it has progesterone receptors. A smallerpercentage of breast cancers, about 20%, may have an excess amount of HER2protein (HER2+). This means that the cancer cells may receive signals from ER,PR or HER2 to promote their growth. Hormonal therapies can disrupt the activ-ity of these hormones to slow down the growth rate or even stop the proliferation.Triple-negative breast cancer (TNBC) is the one that gives negative test result forestrogen receptors, progesterone receptors, and excess HER2 protein. Because ofthe lack of these hormones, it is likely that the patients with TNBC do not respondto hormone therapies.Tumour samples from a HER2+ patient, and a Triple-negative breast cancerpatient were taken and transplanted in mouse to generate patient-deriven xenografts(PDX) across multiple generations (passages). Sample SA532 is from a HER2+patient, and it was serially passaged for up to 10 generations, X1, X2, . . ., X10.The sequenced data compromises 8381 cells. Sample SA609 is from a patient withTriple-negative breast cancer, and it was serially passaged for up to 10 generations,X1, X2, . . ., X10. The sequenced data compromises the total number of 10553cells. We have used the cells from the passages that is studied by Salehi et al.834.3 Data analysis4.3.1 Data cleaningThe initial cleaning steps remove cells with the number of read counts less than500,000 or the quality below 0.75 [76]. Also the s-phase cells are removed at thisstep as described by Salehi et al. [104]. The total number of 1345 cells are retainedfrom OV 2295. However, we only used 891 cells from the study reported by Lakset al. in [76]. As a first approximation to clones, Laks et al. used dimensionalityreduction and clustering to identify 9 cell subsets with shared copy number profiles.Then 891 cells are selected from subsets with more than 50 cells. The number ofcells retained after the initial cleaning for SA532 and SA609 is 2400 and 3243,respectively.4.3.2 Main analysisThe first step in the analysis is to prepare the copy number data and the underlyingtree using the aligned reads. A schematic view of input data is depicted in Fig-ure 4.1(a). The copy number of the cells is inferred using the approach suggestedby Zahn et. al in [141]. As a result, the whole genome of each cell is segmentedby bins, and each bin is associated with a copy number and a quality measure. Aschematic view of the CNV data is depicted in Figure 4.1(b). The bins with qualityless than 0.99 are ignored in the Corrupt analysis. Corrupt utilizes the CNV callbins to infer the underlying tree (See Figure 4.1 (c)) [1]. Corrupt is then run with10,000 number of iterations across 10 chains. The false positive and false nega-tive error rates are set globally, and the maximum rate for false positive and falsenegative are 0.1 and 0.5, respectively [1].Figures 4.5, 4.6 and 4.7 depict the copy number profile and the inferred treeusing Corrupt. The tree is cut to identify subset of cells sharing similar copy num-ber profiles by using the approach of Salehi et al. in [104]. For SA609, we usedthe super tree inferred by Salehi et al., and removed the cells that are not in the setof SA609 cells. The distribution of the number of cells from different samples ineach clone for OV 2295, SA532, and SA609 is depicted in Figures 4.2, 4.3 and 4.4,respectively.84Cell 1 Cell 2 Cell 3 Cell M Cell i … … Cell 4 (a) CNV data Bin 1 Bin 2 Bin 3 Bin 4 Bin 5 … Bin j-1 Bin j Bin j+1 … Bin Z Cell 1 … Cell i ci,j … Cell M (b) _Locus 1 Locus 2 Locus 3 Locus 4 … Locus i … Locus N (d) Read Counts Locus 1 Locus 2 Locus 3 Locus 4 … Locus j … Locus N Cell 1 … Cell i (bi,j, di,j) … Cell M (e) SNV probability Locus 1 Locus 2 Locus 3 Locus 4 … Locus j … Locus N Cell 1 … Cell i pi,j … Cell M (f) (c) (a) (b) (c) (e) (d) (f) CNV inference SNV inference Pseudo-bulk SNV detection Counts data Tree inference (g) Figure 4.1: (Caption next page.)85Figure 4.1: (a) A schematic view of single cell sequencing data. The blackshort horizontal lines represent reads that are associated with each cell’sidentity. The vertical grey lines separate different segments of thegenome. (b) The CNV data inferred by using the single cell sequenc-ing data depicted in (a). The genome is segmented into continuous bins.The CNV calling algorithm assigns copy number ci, j to cell i and binj. (c) A hypothetical phylogenetic tree inferred by Corrupt. The blackcircles represent cells, and the blank circles represent CNV alteration.(d) The list of loci with high SNV call probability from the pseduo-bulkanalysis. (e) The raw data, (bi, j,di, j), denote the number of variant andmatched reads at loci j and cell i. (f) SNV call probability inferred byCellMutScope, for loci j and cell i, is denoted by pi, j. (g) The depen-dencies between different parts of the analysis.A B C D E F G HCNV tree clusters0.00.10.20.30.40.50.6Cells precentageOV2295OV2295(R)TOV2295(R)Figure 4.2: The distribution of the number of cells from different samples ineach clone for OV 2295.86A B C DCNV tree clusters0.00.20.40.60.81.0Cells precentageX1X2X3X4X5X7X8X9Figure 4.3: The distribution of the number of cells from different samples ineach clone for SA532.A B C D E FCNV tree clusters0.00.10.20.30.40.50.60.70.8Cells precentageX1X2X3X4X5X6X7X8X9Figure 4.4: The distribution of the number of cells from different samples ineach clone for SA609.871 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 X YCloneDCABEFIHGSampleCloneABCDEFGHISampleSA921SA922SA1090Copy Number01234567891011+Figure 4.5: The copy number heatmap for OV 2295 data. The inferred underlying phylogeny tree based on CNV data ison the left side of the heatmap. The samples and the clones are colour coded. The cells are from three samplesOV 2295 (SA1090), OV 2295(R) (SA921) and TOV 2295(R) (SA922). The tree includes 9 clones identified in Lakset al. study [76].881 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 X YCloneAUnBDCSampleCloneABCDUnSampleSA532X1XB00118SA532X2XB00147SA532X3XB00210SA532X4XB00273SA532X5XB00478SA532X7XB01216SA532X8XB01398SA532X9XB02755Copy Number01234567891011+Figure 4.6: The copy number heatmap for SA532 data. The inferred underlying phylogeny tree based on CNV data ison the left side of the heatmap. The samples and the clones are colour coded. The clones are identified using theapproach by Salehi et al. study [104].891 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 X YCloneBACDUnEFGHSampleCloneABCDEFGHUnSampleSA609X1XB00290SA609X2XB00426SA609X3XB01584SA609X4XB01721SA609X5XB01844SA609X6XB01899SA609X7XB02184SA609X8XB02312SA609X9XB02360SA609X10XB02454Copy Number01234567891011+Figure 4.7: The copy number heatmap for SA609 data. The inferred underlying phylogeny tree based on CNV data ison the left side of the heatmap. The samples and the clones are colour coded. The clones are identified using theapproach by Salehi et al. study [104].90Next, a group of SNV loci candidates is selected from the aggregated data inwhich all cells were collapsed into a ‘pseudo-bulk’ genome (See Figure 4.1 (d)).MutationSeq and Strelka are run for selecting SNV loci candidates from the pseudo-bulk sample [31, 105]. The loci with MutationSeq probability greater than 0.9and Strelka’s Phred quality score of 20 or more are identified as SNV candidates.In addition, each locus should demonstrate at least two variant reads across thecells. The total number of 14020, 11416 and 15446 SNV loci are selected forOV 2295, SA532 and SA609, respectively. The cells with no variant across the lociare removed. This condition reduces the number of cells to 2400 for SA532, 3243for SA609, and 731 for OV 2295. Table 4.1 summarizes some statistics about thecells used in the following analysis. As specified in this table, OV 2295 has a largernumber of reads, and a higher breadth, and a higher depth coverage in comparisonto SA532 and SA609.CellMutScope described in Chapter 3, incorporates the reads data (See Fig-ure 4.1 (e)), the copy number data (See Figure 4.1 (b)), and the underlying tree(See Figure 4.1 (c)) to detect SNVs per cell. Output of CellMutScope is schemati-cally depicted in Figure 4.1 (g). The value of εFN is approximated asεFN =ndeletionntotal× 12, (4.1)where ndeletion is the number of bins with copy number less than 2, and ntotal isthe total number of bins across the genome. The value 12 is the probability thatthe mutation has occurred on the missing strand. The value of εFP is 0.001. Theoverall picture of the pipeline and dependencies is depicted in Figure 4.1 (g).4.4 Experimental resultsCellMutScope is applied on OV 2295, SA532, and SA609. Since the model is basedon placing SNV loci on the tree, the distribution of SNV call probabilities is bi-modal. It depends on the position of the node (if it is an ancestor of a cell or not).Figures 4.8, 4.9 and 4.10 depict the SNV call probability distributions for OV 2295,SA532 and SA609, respectively. Figures 4.8 (a), 4.9 (a) and 4.10 (a) confirm thatthe distribution has two peaks at 0 and 1. Zoomed in at the distribution of SNV call91(b)(a)Figure 4.8: (a) SNV call probability distribution for OV 2295. (b) SNV callprobability distribution for OV 2295 for the probabilities between 0.05and 0.95.probabilities between 0.05 and 0.95 are depicted in Figures 4.8 (b), 4.9 (b) and 4.10(b). There is also a peak at 0.5 that is smaller than the peaks at 0 and 1. The 0.5peak represents for the loci that the model was not confident in SNV calling. Thiscan be partly explained by missing data, low coverage, or inconsistency of SNVdata with the tree. The coverage data supports this assumption as well.Figure 4.11 depicts the SNV call probabilities for OV 2295. For the sake ofvisualization, a list of 200 SNV loci candidates is randomly sampled, from the setof all loci with equal probabilities. Moreover, the loci are sorted using R’s ‘hclust’function with ‘average’ metric and ‘manhattan’ distance. Figure 4.12 depicts thenumber of variant reads in each cell across the same set of loci. The loci in SNVcall probability plot (Figure 4.11) and the plot that depicts the number of variantreads (Figure 4.12) follow a exactly similar order. Considering both Figures 4.12and 4.11, the number of variant reads in the area (black area in Figure 4.12) thatCellMutScope has low probability (blue area in Figure 4.11) is small.92(a) (b)Figure 4.9: (a) SNV call probability distribution for SA532. (b) SNV call prob-ability distribution for SA532 for the probabilities between 0.05 and0.95.(a) (b)Figure 4.10: (a) SNV call probability distribution for SA609. (b) SNV callprobability distribution for SA609 for the probabilities between 0.05and 0.95.93Metric OV2295 SA609 SA532Number of cells 731 3243 2400Number of SNV loci candidates 14020 15446 11416Total reads 8.558316e+06 1.446370e+06 2.635680e+06Coverage breadth 1.563200e-01 3.452980e-02 5.356769e-02Coverage depth 1.776126e-01 3.615207e-02 5.871972e-02Total mapped reads 7.913032e+06 1.329028e+06 2.376801e+06Unmapped reads 6.452840e+05 1.173421e+05 2.588786e+05Table 4.1: The table shows the total number of cells, the total number of SNVloci candidates, the average number of reads, the average breadth cov-erage, the average depth coverage, the average total number of mappedreads and the average total number of unmapped reads across the cells ineach sample.94CloneDCABEFIHGSampleCloneABCDEFGHISampleSA921SA922SA1090SNV call probabilities00.20.40.60.81Figure 4.11: The heatmap depicts OV 2295 SNV call probability of 200 random loci across the cells. The blue colourdepicts the low probabilities near 0, the white colour depicts the probabilities around 0.5, and the orange colourdepicts the probabilities around 1. The inferred underlying phylogeny tree based on CNV data is on the leftside of the heatmap. The samples and the clones are colour coded. The cells are from three samples OV 2295(SA1090), OV 2295(R) (SA921) and TOV 2295(R) (SA922). The tree includes 9 clones identified in Laks et al.study [76].95CloneDCABEFIHGSampleCloneABCDEFGHISampleSA921SA922SA1090Number of variant allele012345Figure 4.12: The heatmap depicts the number of variant reads observed in each cell of OV 2295 for 200 random loci.No variant read is depicted in black, and reads with variant alleles are depicted in bright colours. The inferredunderlying phylogeny tree based on CNV data is on the left side of the heatmap. The samples and the clonesare colour coded. The cells are from three samples OV 2295 (SA1090), OV 2295(R) (SA921) and TOV 2295(R)(SA922). The tree includes 9 clones identified in Laks et al. study [76].96For SA532, Figure 4.13 depicts the SNV call probabilities, and Figure 4.14depicts the number of variant reads for 200 selected loci across the cells. A similarapproach is used for sorting SNV loci across the x-axis. As depicted in Figure 4.13,more detailed clonal structure appears in clone A that was not apparent with onlythe CNV data and smaller group of cells sharing similar SNV profile. Comparedwith the results from OV 2295, we see more variant reads in the area (black areain Figure 4.14) that CellMutScope has low probability (blue area in Figure 4.13).For SA609, Figures 4.15 and 4.16 show SNV call probabilities and the number ofvariant read across 200 random loci. As shown, subgroups of cells in clone Cand H sharing similar SNVs can represent sub-clones in clone C and H. Similarly,for SA609 (see Figures 4.15 and 4.16), there are a few cells having variant readswhile they do not have a high SNV call probability. This suggests the need for aquality control step to remove the inconsistent loci (i.e. the loci having SNV datainconsistent with the tree structure). Comparing the results from OV 2295, SA532and SA609, we can see SA532 and SA609 have a higher rate of cells with variantreads while they do not have a high SNV call probability. This can be explainedby the difference between the coverage depth of the samples. OV 2295 has highercoverage depth. The coverage depth of OV 2295, SA532 and SA609 is 0.17, 0.03and 0.05, respectively (See Table 4.1).Quality control step A SNV is called present, if the value of SNV call probabilityfrom CellMutScope is greater or equals to 0.8. For clones that SNV is not present,we count the number of cells with variant reads normalized to the total number ofcells in the clone. The average number of variant reads is expected to be greaterthan the average coverage times the number of cells in each clone divided by thecopy number (assuming only one of the copies has mutation). This implies thenormalized number of variant reads should be greater than the average coveragedivided by copy number. With the average coverage of 0.05 and 0.01 and copynumber 10, the threshold is 0.005 and 0.001, respectively.97CloneAUnBDCSampleCloneABCDUnSampleSA532X1XB00118SA532X2XB00147SA532X3XB00210SA532X4XB00273SA532X5XB00478SA532X7XB01216SA532X8XB01398SA532X9XB02755SNV call probabilities00.20.40.60.81Figure 4.13: The heatmap depicts SA532 SNV call probability for 200 random loci across the cells. The blue colourdepicts the low probabilities near 0, the white colour depicts the probabilities around 0.5, and the orange colourdepicts the probabilities around 1. The inferred underlying phylogeny tree based on CNV data is on the left sideof the heatmap. The samples and the clones are colour coded. The clones are identified using the approach bySalehi et al. study [104].98CloneAUnBDCSampleCloneABCDUnSampleSA532X1XB00118SA532X2XB00147SA532X3XB00210SA532X4XB00273SA532X5XB00478SA532X7XB01216SA532X8XB01398SA532X9XB02755Number of variant allele01234567810Figure 4.14: The heatmap depicts the number of variant reads observed in each cell of SA532 for 200 random loci.No variant read is depicted in black, and reads with variant alleles are depicted in bright colours. The inferredunderlying phylogeny tree based on CNV data is on the left side of the heatmap. The samples and the clones arecolour coded. The clones are identified using the approach by Salehi et al. study [104].99CloneBACDUnEFGHSampleCloneABCDEFGHUnSampleSA609X1XB00290SA609X2XB00426SA609X3XB01584SA609X4XB01721SA609X5XB01844SA609X6XB01899SA609X7XB02184SA609X8XB02312SA609X9XB02360SA609X10XB02454SNV call probabilities00.20.40.60.81Figure 4.15: The heatmap depicts SA609 SNV call probability for 200 random loci across the cells. The blue colourdepicts the low probabilities near 0, the white colour depicts the probabilities around 0.5, and the orange colourdepicts the probabilities around 1. The inferred underlying phylogeny tree based on CNV data is on the left sideof the heatmap. The samples and the clones are colour coded. The clones are identified using the approach bySalehi et al. study [104].100CloneBACDUnEFGHSampleCloneABCDEFGHUnSampleSA609X1XB00290SA609X2XB00426SA609X3XB01584SA609X4XB01721SA609X5XB01844SA609X6XB01899SA609X7XB02184SA609X8XB02312SA609X9XB02360SA609X10XB02454Number of variant allele0123456810Figure 4.16: The heatmap depicts the number of variant reads observed in each cell of SA609 for 200 random loci.No variant read is depicted in black, and reads with variant alleles are depicted in bright colours. The inferredunderlying phylogeny tree based on CNV data is on the left side of the heatmap. The samples and the clones arecolour coded. The clones are identified using the approach by Salehi et al. study [104].101CloneDCABEFIHGSampleCloneABCDEFGHISampleSA921SA922SA1090SNV call probabilities00.20.40.60.81Figure 4.17: The heatmap depicts OV 2295 SNV call probability for 200 random loci passed the quality control step withthreshold equals 0.001 across the cells. The blue colour depicts the low probabilities near 0, the white colourdepicts the probabilities around 0.5, and the orange colour depicts the probabilities around 1. The inferredunderlying phylogeny tree based on CNV data is on the left side of the heatmap. The samples and the clonesare colour coded. The cells are from three samples OV 2295 (SA1090), OV 2295(R) (SA921) and TOV 2295(R)(SA922). The tree includes 9 clones identified in Laks et al. study [76].102CloneDCABEFIHGSampleCloneABCDEFGHISampleSA921SA922SA1090Number of variant allele012345Figure 4.18: The heatmap depicts the number of variant reads observed in each cell of OV 2295 for 200 random locipassed the quality control step with threshold equals 0.001 across the cells. No variant read is depicted in black,and reads with variant alleles are depicted in bright colours. The inferred underlying phylogeny tree based onCNV data is on the left side of the heatmap. The samples and the clones are colour coded. The cells are fromthree samples OV 2295 (SA1090), OV 2295(R) (SA921) and TOV 2295(R) (SA922). The tree includes 9 clonesidentified in Laks et al. study [76].103CloneDCABEFIHGSampleCloneABCDEFGHISampleSA921SA922SA1090SNV call probabilities00.20.40.60.81Figure 4.19: The heatmap depicts the SNV call probability of 200 random loci passed the quality control step withthreshold equals 0.005 across the cells. The blue colour depicts the low probabilities near 0, the white colourdepicts the probabilities around 0.5, and the orange colour depicts the probabilities around 1. The inferredunderlying phylogeny tree based on CNV data is on the left side of the heatmap. The samples and the clonesare colour coded. The cells are from three samples OV 2295 (SA1090), OV 2295(R) (SA921) and TOV 2295(R)(SA922). The tree includes 9 clones identified in Laks et al. study [76].104CloneDCABEFIHGSampleCloneABCDEFGHISampleSA921SA922SA1090Number of variant allele012345Figure 4.20: The heatmap depicts the number of variant reads observed in each cell for 200 random loci passed thequality control step with threshold equals 0.005 across the cells. No variant read is depicted in black, and readswith variant alleles are depicted in bright colours. The inferred underlying phylogeny tree based on CNV data ison the left side of the heatmap. The samples and the clones are colour coded. The cells are from three samplesOV 2295 (SA1090), OV 2295(R) (SA921) and TOV 2295(R) (SA922). The tree includes 9 clones identified inLaks et al. study [76].105Sample Original TH=0.005 TH=0.001OV2295 14020 12574 12574SA609 15446 13052 9962SA532 11416 9648 8403Table 4.2: The table shows the total number of cells before quality controlstep (original) and the total number of cells after quality control withthreshold equal 0.005 and 0.001.For OV 2295, the total number of SNVs passing quality control step with boththresholds 0.005 and 0.001 is 12574 (see Table 4.2). Figures 4.17 and 4.18 depictSNV call probability for 200 loci randomly sampled from the set of loci that pass thequality check with a threshold of 0.001 for OV 2295. Figures 4.19 and 4.20 depictresults for the same data with a threshold of 0.005. As shown in Table 4.2, forSA532, the total number of SNVs passing quality control step with thresholds 0.005and 0.001 are 9648 and 8403, respectively. For SA609, the numbers are 13052and 9962 for thresholds of 0.005 and 0.001, respectively. The results of similaranalysis for SA532 and SA609 with thresholds of 0.001 and 0.005 are depicted inFigures 4.21 to 4.28. See Figures 4.21 and 4.22 for SA532 results with a thresholdof 0.001. See Figures 4.23 and 4.24 for SA532 results with a threshold of 0.005.See Figures 4.25 and 4.26 for SA609 results with a threshold of 0.001. See Figures4.27 and 4.28 for SA609 results with a threshold of 0.005. As shown in the resultsafter the quality control step, we observe lower rate of cells having variant readswhile they do not have a high SNV call probability.106CloneAUnBDCSampleCloneABCDUnSampleSA532X1XB00118SA532X2XB00147SA532X3XB00210SA532X4XB00273SA532X5XB00478SA532X7XB01216SA532X8XB01398SA532X9XB02755SNV call probabilities00.20.40.60.81Figure 4.21: The heatmap depicts SA532 SNV call probability of 200 random loci passed the quality control step withthreshold equals 0.001. The blue colour depicts the low probabilities near 0, the white colour depicts the proba-bilities around 0.5, and the orange colour depicts the probabilities around 1. The inferred underlying phylogenytree based on CNV data is on the left side of the heatmap. The samples and the clones are colour coded. Theclones are identified using the approach by Salehi et al. study [104].107CloneAUnBDCSampleCloneABCDUnSampleSA532X1XB00118SA532X2XB00147SA532X3XB00210SA532X4XB00273SA532X5XB00478SA532X7XB01216SA532X8XB01398SA532X9XB02755Number of variant allele01234567810Figure 4.22: The heatmap depicts the number of variant reads observed in each cell of SA532 for 200 random locipassed the quality control step with threshold equals 0.001. No variant read is depicted in black, and reads withvariant alleles are depicted in bright colours. The inferred underlying phylogeny tree based on CNV data is onthe left side of the heatmap. The samples and the clones are colour coded. The clones are identified using theapproach by Salehi et al. study [104].108CloneAUnBDCSampleCloneABCDUnSampleSA532X1XB00118SA532X2XB00147SA532X3XB00210SA532X4XB00273SA532X5XB00478SA532X7XB01216SA532X8XB01398SA532X9XB02755SNV call probabilities00.20.40.60.81Figure 4.23: The heatmap depicts SA532 SNV call probability of 200 random loci passed the quality control step withthreshold equals 0.005. The blue colour depicts the low probabilities near 0, the white colour depicts the proba-bilities around 0.5, and the orange colour depicts the probabilities around 1. The inferred underlying phylogenytree based on CNV data is on the left side of the heatmap. The samples and the clones are colour coded. Theclones are identified using the approach by Salehi et. al. study [104].109CloneAUnBDCSampleCloneABCDUnSampleSA532X1XB00118SA532X2XB00147SA532X3XB00210SA532X4XB00273SA532X5XB00478SA532X7XB01216SA532X8XB01398SA532X9XB02755Number of variant allele0123456791011+Figure 4.24: The heatmap depicts the number of variant reads observed in each cell of SA532 for 200 random locipassed the quality control step with threshold equals 0.005. No variant read is depicted in black, and reads withvariant alleles are depicted in bright colours. The inferred underlying phylogeny tree based on CNV data is onthe left side of the heatmap. The samples and the clones are colour coded. The clones are identified using theapproach by Salehi et al. study [104].110CloneBACDUnEFGHSampleCloneABCDEFGHUnSampleSA609X1XB00290SA609X2XB00426SA609X3XB01584SA609X4XB01721SA609X5XB01844SA609X6XB01899SA609X7XB02184SA609X8XB02312SA609X9XB02360SA609X10XB02454SNV call probabilities00.20.40.60.81Figure 4.25: The heatmap depicts SA609 SNV call probability of 200 random loci passed the quality control step withthreshold equals 0.001. The blue colour depicts the low probabilities near 0, the white colour depicts the proba-bilities around 0.5, and the orange colour depicts the probabilities around 1. The inferred underlying phylogenytree based on CNV data is on the left side of the heatmap. The samples and the clones are colour coded. Theclones are identified using the approach by Salehi et al. study [104].111CloneBACDUnEFGHSampleCloneABCDEFGHUnSampleSA609X1XB00290SA609X2XB00426SA609X3XB01584SA609X4XB01721SA609X5XB01844SA609X6XB01899SA609X7XB02184SA609X8XB02312SA609X9XB02360SA609X10XB02454Number of variant allele01234568910Figure 4.26: The heatmap depicts the number of variant reads observed in each cell of SA609 for 200 random locipassed the quality control step with threshold equals 0.001. No variant read is depicted in black, and reads withvariant alleles are depicted in bright colours. The inferred underlying phylogeny tree based on CNV data is onthe left side of the heatmap. The samples and the clones are colour coded. The clones are identified using theapproach by Salehi et al. studies [104].112CloneBACDUnEFGHSampleCloneABCDEFGHUnSampleSA609X1XB00290SA609X2XB00426SA609X3XB01584SA609X4XB01721SA609X5XB01844SA609X6XB01899SA609X7XB02184SA609X8XB02312SA609X9XB02360SA609X10XB02454SNV call probabilities00.20.40.60.81Figure 4.27: The heatmap depicts SA609 SNV call probability of 200 random loci passed the quality control step withthreshold equals 0.005. The blue colour depicts the low probabilities near 0, the white colour depicts the proba-bilities around 0.5, and the orange colour depicts the probabilities around 1. The inferred underlying phylogenytree based on CNV data is on the left side of the heatmap. The samples and the clones are colour coded. Theclones are identified using the approach by Salehi et al. study [104].113CloneBACDUnEFGHSampleCloneABCDEFGHUnSampleSA609X1XB00290SA609X2XB00426SA609X3XB01584SA609X4XB01721SA609X5XB01844SA609X6XB01899SA609X7XB02184SA609X8XB02312SA609X9XB02360SA609X10XB02454Number of variant allele012345678910Figure 4.28: The heatmap depicts the number of variant reads observed in each cell of SA609 for 200 random locipassed the quality control step with threshold equals 0.005. No variant read is depicted in black, and reads withvariant alleles are depicted in bright colours. The inferred underlying phylogeny tree based on CNV data is onthe left side of the heatmap. The samples and the clones are colour coded. The clones are identified using theapproach by Salehi et al. study [104].114To confirm that the structure of the tree is an important piece of informationfor detecting the SNVs, particularly the ones that are only present in a subset ofcells (clone specific), we introduce a perturbed tree by permuting the cells on thetree while the tree structure is fixed. We repeated this process 20 times. Fig-ures 4.29, 4.31 and 4.33 show the mean probability of SNV calls across the per-turbed trees. The results show that the perturbation disables the model from de-tecting the clone specific mutations, and only the ancestral mutations survived.Figures 4.30, 4.32 and 4.34 show the SNV variant reads corresponding to the per-mutation analysis of Figures 4.29, 4.31 and 4.33, respectively. This proves the treestructure is a crucial element of the model to enable capturing the clone specificmutations and addressing the sparsity of the data.115CloneDCABEFIHGSampleCloneABCDEFGHISampleSA921SA922SA1090SNV call probabilities00.20.40.60.81Figure 4.29: The heatmap depicts the mean SNV call probability of 200 random loci across the cells of OV 2295 over20 perturbed tree (quality control step with threshold equals 0.001 is applied). The blue colour depicts the lowprobabilities near 0, the white colour depicts the probabilities around 0.5, and the orange colour depicts theprobabilities around 1. The inferred underlying phylogeny tree based on CNV data is on the left side of theheatmap. The samples and the clones are colour coded. The cells are from three samples OV 2295 (SA1090),OV 2295(R) (SA921) and TOV 2295(R) (SA922). The tree includes 9 clones identified in Laks et al. study [76].116CloneDCABEFIHGSampleCloneABCDEFGHISampleSA921SA922SA1090Number of variant allele0123456Figure 4.30: The heatmap depicts the number of variant reads observed in each cell of OV 2295 for 200 random locipassed the quality control step with threshold equals 0.001. No variant read is depicted in black, and reads withvariant alleles are depicted in bright colours. The inferred underlying phylogeny tree based on CNV data is on theleft side of the heatmap. The samples and the clones are colour coded. The cells are from three samples OV 2295(SA1090), OV 2295(R) (SA921) and TOV 2295(R) (SA922). The tree includes 9 clones identified in Laks et al.study [76].117CloneAUnBDCSampleCloneABCDUnSampleSA532X1XB00118SA532X2XB00147SA532X3XB00210SA532X4XB00273SA532X5XB00478SA532X7XB01216SA532X8XB01398SA532X9XB02755SNV call probabilities00.20.40.60.81Figure 4.31: The heatmap depicts the mean SNV call probability of 200 random loci across the cells of SA532 over20 perturbed tree (quality control step with threshold equals 0.001 is applied). The blue colour depicts the lowprobabilities near 0, the white colour depicts the probabilities around 0.5, and the orange colour depicts theprobabilities around 1. The inferred underlying phylogeny tree based on CNV data is on the left side of theheatmap. The samples and the clones are colour coded. The clones are identified using the approach by Salehiet al. study [104].118CloneAUnBDCSampleCloneABCDUnSampleSA532X1XB00118SA532X2XB00147SA532X3XB00210SA532X4XB00273SA532X5XB00478SA532X7XB01216SA532X8XB01398SA532X9XB02755Number of variant allele01234567810Figure 4.32: The heatmap depicts the number of variant reads observed in each cell of SA532 for 200 random locipassed the quality control step with threshold equals 0.001. No variant read is depicted in black, and reads withvariant alleles are depicted in bright colours. The inferred underlying phylogeny tree based on CNV data is onthe left side of the heatmap. The samples and the clones are colour coded. The clones are identified using theapproach by Salehi et al. studies [104].119CloneBACDUnEFGHSampleCloneABCDEFGHUnSampleSA609X1XB00290SA609X2XB00426SA609X3XB01584SA609X4XB01721SA609X5XB01844SA609X6XB01899SA609X7XB02184SA609X8XB02312SA609X9XB02360SA609X10XB02454SNV call probabilities00.20.40.60.81Figure 4.33: The heatmap depicts the mean SNV call probability of 200 random loci across the cells of SA609 over20 perturbed tree (quality control step with threshold equals 0.001 is applied). The blue colour depicts the lowprobabilities near 0, the white colour depicts the probabilities around 0.5, and the orange colour depicts theprobabilities around 1. The inferred underlying phylogeny tree based on CNV data is on the left side of theheatmap. The samples and the clones are colour coded. The clones are identified using the approach by Salehiet al. study [104].120CloneBACDUnEFGHSampleCloneABCDEFGHUnSampleSA609X1XB00290SA609X2XB00426SA609X3XB01584SA609X4XB01721SA609X5XB01844SA609X6XB01899SA609X7XB02184SA609X8XB02312SA609X9XB02360SA609X10XB02454Number of variant allele0123456810Figure 4.34: The heatmap depicts the number of variant reads observed in each cell of SA609 for 200 random locipassed the quality control step with threshold equals 0.001. No variant read is depicted in black, and reads withvariant alleles are depicted in bright colours. The inferred underlying phylogeny tree based on CNV data is onthe left side of the heatmap. The samples and the clones are colour coded. The clones are identified using theapproach by Salehi et al. study [104].1214.5 Mutations in high impact genesWe first annotate (using SnpEff [23]) the SNV loci that pass the quality controlstep. The loci that are annotated as a high impact gene (or a non-synonymouscoding genes) and are only present in at most two clones are identified. The clonespecific genes for OV 2295, SA532 and SA609 are listed in Tables 4.3, 4.4 and 4.5,respectively. For OV 2295, there are SNVs in GLG1 that are specific to clone C,and there are SNVs in DNHD1 that are specific to clone H. For SA532, there areSNVs in BCCIP, PTDSS2, PLEKHH1, UBALD1, NEURL4, MGAT5B, DOCK6,CBLC, TRIM28, PARS2, TP53RK, FKBP7, SPEG, TRANK1, SEL1L3, GRIA1,HIST1H2AM, PSMB9 and DHRS13 genes that are only present in clone A. SomeSNVs that are only present in clone B occur on SMARCA5; and some of the onesthat are only present in clone C occur on NOC3L and SPEG. Moreover, SNVs onXIAP gene are only present in clone D. For SA609, there are SNVs that are onlypresent in clone D and occur in TLX1, OTUD7A, ZNF112 and MROH2A; andthere are SNVs only present in clone E and occur in ARHGAP21, SLC12A5 andZNF622. Some of the clone H specific mutations occur in DGKA, CARHSP1,PANK4, TTN, LRRC61, BAI1, MAGEC3 and NHS.We explore if any of these genes are in cancer genes census. We find TLX1 inclone D of SA609 and CBLC in clone A of SA532 in cancer gene census. Also,we find SMARCA5 in clone B of SA532, and SLC34A2 in clone B and C ofSA609. Both of the genes are reported as translocation partner genes. Reviewingthe literature, it is interesting that we find studies reporting the above genes effec-tive in different cancers. Hideshima et al. report TP53RK confers poor prognosisin multiple myeloma tumours. Huang et al. report XIAP possesses a critical rolein promotion of cell survival and maintenance of cellular homeostasis for breastand colon carcinoma [57]. BCCIP is reported to have a role in the maintenanceof genomic integrity [58, 85]. Jin et al. report overexpression of SMARCA5 cor-relates with cell proliferation and migration in breast cancer [62]. Yu et al. reportSLC12A5 as an oncogene in clone cancer [138]. Further gene expression analysiscan be done to confirm if any of the identified mutations has a significant impact inthe dynamics of studied tumours.122Clone Gene listC GLG1H DNHD1Table 4.3: High impact and non-synonymous coding genes identified havingSNVs that are only present in at most 2 clones of OV 2295.Clone Gene listA BCCIP, PTDSS2, PLEKHH1, UBALD1, NEURL4, GPX4GRIA1, HIST1H2AM, PSMB9, DHRS13, HSPB6MGAT5B, DOCK6, CBLC, TRIM28, PARS2,TP53RK, FKBP7, SPEG, TRANK1, SEL1L3,B SMARCA5C NOC3L, SPEGD XIAPB, D PYCRL , CCKBR, SFI1, ANK2Table 4.4: High impact and non-synonymous coding genes identified havingSNVs that are only present in at most 2 clones of SA532.4.6 ConclusionWe applied CellMutScope to three real datasets each compromising thousands ofclonally related cells. The model detects SNVs per cell from the input data withvery low breadth and depth coverage (See Table 4.2). Although the results showthat better performance is achieved with a higher coverage, the applicable coverageis still lower than the ones achieved in other existing methods. In order to havea DLP data with higher coverage we can either increase the number of lanes orsynthetically group similar cells as a pseudo-cell with higher coverage. The latterapproach can introduce some errors as cells with different genomic profiles maybemistakenly grouped. The results demonstrates the CellMutScope identifies clonespecific mutations that are typically considered hardly detectable. In some of theexisting clones, minor clonal structures (that were not distinguished with only theCNV data) are identified. A few of clone specific SNVs occur in high impact andnon-synonymous coding genes. This may suggest the underling genetic pathways.The next step is to validate if these mutations have real impacts on the expansion123Clone Gene listD TLX1, OTUD7A, ZNF112, MROH2AE ARHGAP21, SLC12A5, ZNF622C, B SHISA2, NOX5, TNFAIP8L1, OTUD7B, MIEF1, CLNK,SLC34A2E, G OTUD1, GCGR, MON1A, RASSF1,CWH43, DNAJC25, CTNNBIP1Table 4.5: High impact and non-synonymous coding genes identified havingSNVs that are only present in at most 2 clones of SA609.or shrinkage of their associated clones. Single cell transcriptome sequencing datashould be exploited to investigate if any of these mutations has any impact on theexpression of the genes.124Chapter 5Conclusion“Learn from yesterday, live for today, hope for tomorrow. Theimportant thing is not to stop questioning.”- Albert Einstein5.1 Summary of contributionsThis dissertation outlines two computational methods and their application on realworld cancer datasets. The methods can be used to identify the portfolio of SNVsthrough incorporation of clonal information. The accurate identification of SNVsbenefits the inference of tumour dynamics.Mutation detection and classification through probabilistic integrationof clonal population structureWe developed a statistical method, MuClone, to detect and classify mutations frommultiple clonally related tumour whole genome or exome sequencing samples. Ourmethod models SNV genotype, normal contamination, and clonal prevalence, all ofwhich confound the detection of low prevalence SNVs. We showed that incor-porating clonal information improves the detection of SNVs, particularly the lowprevalence ones. In addition, we confirmed multi sample bulk sequencing data isa viable experimental design (e.g rapid autopsy program that at the time of a pa-125tient’s death, tens to hundreds of metastatic samples are collected for future study)which can be exploited to interpret evolutionary properties of cancer.Single cell somatic mutation detection through incorporation of theunderlying phylogenyWe developed a statistical model, CellMutScope, to detect mutations across thou-sands of clonally related sequenced single cells. Incorporating the underlying phy-logeny and copy number profile of cells, enabled us to overcome the low breadthand depth coverage of the data. By sharing statistical strength, we showed that themodel imputes the missing data and detects SNVs in each cell. We also frame howwe can incorporate SNV data into the Corrupt model (i) to infer the underlying treephylogeny using SNV data with or without CNV data, or (ii) to extend the under-lying CNV data with SNV loci placed on the tree. Therefore, we confirmed thatwith single cell sequencing data, it is possible to reveal the genomic profile of theunderlying tumour clones in more detail.5.2 Future work and discussionThis thesis focused on identification and interpretation of SNVs across clonallyrelated tumour samples while exploiting clonal information. The first proposedmodel in Chapter 2 is restricted to the use of information from a flat clusteringof clones. We believe exploiting all information about the underlying phylogenyof clones will improve the performance of detecting mutations. Although tumourevolution is the result of accumulation of SNVs with other types of genomic aber-rations, the proposed method is only focused on the detection of SNVs. An impor-tant subject for future research is to integrate analysis of different aberrations intumours with the underlying phylogeny to increase accuracy of detecting genomicaberrations, and to decipher tumour evolutionary process. We expect the joint anal-ysis of genomic data will improve the accuracy of detecting genomic aberrations.Exploiting bulk genome sequencing data combined with other type of data (e.g.single cell sequencing data, or transcriptome sequencing data) can provide a morecomprehensive insight.With access to single cell sequencing data, we proposed a model in Chapter 3126to detect SNVs in every individual cell by incorporating the underlying phyloge-netic tree. The model helps distinguish clone specific mutations and it disclosesminor clone structure in the data. Identification of driver SNVs and their fitness canbe extended to include all genomic aberrations. This should enhance the capac-ity of interpretations in tracking various aberrations. In addition, measuring geneexpression of tumour clones links functional consequences to somatic aberrations.Therefore, genomic and transcriptomic single cell sequencing data can comple-ment each other to give us insights on changes of gene expressions as a tumourevolve. This enables us to reconstruct tumours life histories.It is observed that cells evolve through accumulation of genomic aberrations,epigenetic modifications and translation. This results in the heterogeneity of cellpopulations. Single cell sequencing data together with other biological informa-tion, open an avenue for probing the tumour dynamics in various patients by iden-tifying distinct patterns. The difference in tumour dynamic patterns can lead to-wards a novel classification of tumour subtypes. This can help clinical decisionsand treatments. The emerging field is called personalized medicine in which thetreatment is personalized based on an individual’s genomic profile.Our understanding about cancer as a very complex disease can be improved bythe study of the evolutionary process of cells under different conditions. Extensionsof the experimental techniques like cell lines, patient derived xenografted samples,and genetically engineered mouse models can provide data to systematically studythe evolutionary dynamics of cells. For example, we can study the evolutionaryconsequences of introducing mutations into the model or investigate the responseof xenograft populations to therapeutic intervention (like gene knockouts). Themodels presented here can be applied at scale to comprehensively analyse SNVsin different experimental designs. The results can help reveal the mechanisms oftumourigenesis, treatment resistance, and metastatic progression in different con-ditions.In spite of advances in sequencing technologies and different experiments, lit-tle is known about how to predict the evolutionary dynamics of cell populations,and how this information should be used in practice. Little findings from differentexperiments will add up, and will eventually enhance our understanding of tumourdynamics. Therefore, another milestone lies in interrelationship of clinical deci-127sions and computational results. Here, there are critical questions that are must beanswered: What can we learn from inferred tumour dynamics? How should weaggregate computational results for clinically valid conclusions? How many com-putational cases are sufficient for clinical conclusions? Are the discovered tumourdynamic patterns reoccurring? If they are, how long the inferred data will remainvalid in an evolving population? How frequent should we expect to find a de novotumour dynamic pattern? How much do the parameters used in the analysis (e.g.number of samples or distribution of longitudinal/anatomical samples) affect theresults? All these questions reflect the ambiguities and uncertainties laid under thisfield; and they should set our expectations about the impact of the field in practice.I hope the advances in computational methods, in particular the ones we de-scribed here help improve our understanding about cancer and eventually help savemore lives.128Bibliography[1] Corrupt. URL https://github.com/UBC-Stat-ML/nowellpack. → pages55, 57, 59, 60, 63, 65, 73, 84[2] URL https://www.breastcancer.org. → page 83[3] Mutect2. URLhttps://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/org_broadinstitute_gatk_tools_walkers_cancer_m2_MuTect2.php. →pages 13, 35[4] K. Anderson, C. Lutz, F. W. van Delft, C. M. Bateman, Y. Guo, S. M.Colman, H. Kempski, A. V. Moorman, I. Titley, J. Swansbury, L. Kearney,T. Enver, and M. Greaves. Genetic variegation of clonal architecture andpropagating cells in leukaemia. Nature, 469(7330):356–61, 2011. → page3[5] L. Angus, M. Smid, S. M. Wilting, J. van Riet, A. Van Hoeck, L. Nguyen,S. Nik-Zainal, T. G. Steenbruggen, V. C. G. Tjan-Heijnen, M. Labots, J. M.G. H. van Riel, H. J. Bloemendal, N. Steeghs, M. P. Lolkema, E. E. Voest,H. J. G. van de Werken, A. Jager, E. Cuppen, S. Sleijfer, and J. W. M.Martens. The genomic landscape of metastatic breast cancer highlightschanges in mutation and signature frequencies. Nature Genetics, 51(10):1450–1458, 2019. → page 82[6] S. A. Aparicio and D. G. Huntsman. Does massively parallel dnaresequencing signify the end of histopathology as we know it? The Journalof Pathology: A Journal of the Pathological Society of Great Britain andIreland, 220(2):307–315, 2010. → page 8[7] C. E. Barbieri, S. C. Baca, M. S. Lawrence, F. Demichelis, M. Blattner,J.-P. Theurillat, T. A. White, P. Stojanov, E. Van Allen, N. Stransky,E. Nickerson, S.-S. Chae, G. Boysen, D. Auclair, R. C. Onofrio, K. Park,129N. Kitabayashi, T. Y. MacDonald, K. Sheikh, T. Vuong, C. Guiducci,K. Cibulskis, A. Sivachenko, S. L. Carter, G. Saksena, D. Voet, W. M.Hussain, A. H. Ramos, W. Winckler, M. C. Redman, K. Ardlie, A. K.Tewari, J. M. Mosquera, N. Rupp, P. J. Wild, H. Moch, C. Morrissey, P. S.Nelson, P. W. Kantoff, S. B. Gabriel, T. R. Golub, M. Meyerson, E. S.Lander, G. Getz, M. A. Rubin, and L. A. Garraway. Exome sequencingidentifies recurrent SPOP, FOXA1 and MED12 mutations in prostatecancer. Nature Genetics, 44(6):685–689, 2012. → page 82[8] A. Bashashati, G. Ha, A. Tone, J. Ding, L. M. Prentice, A. Roth, J. Rosner,K. Shumansky, S. Kalloger, J. Senz, W. Yang, M. McConechy, N. Melnyk,M. Anglesio, M. T. Y. Luk, K. Tse, T. Zeng, R. Moore, Y. Zhao, M. A.Marra, B. Gilks, S. Yip, D. G. Huntsman, J. N. McAlpine, and S. P. Shah.Distinct evolutionary trajectories of primary high-grade serous ovariancancers revealed through spatial mutational profiling. Journal of Pathology,231(1):21–34, 2013. → page 3[9] T. Baslan, J. Kendall, L. Rodgers, H. Cox, M. Riggs, A. Stepansky,J. Troge, K. Ravi, D. Esposito, B. Lakshmi, M. Wigler, N. Navin, andJ. Hicks. Genome-wide copy number analysis of single cells. NatureProtocols, 7(6):1024–1041, 2012. → page 56[10] T. Baslan, J. Kendall, B. Ward, H. Cox, A. Leotta, L. Rodgers, M. Riggs,S. D’Italia, G. Sun, M. Yong, K. Miskimen, H. Gilmore, M. Saborowski,N. Dimitrova, A. Krasnitz, L. Harris, M. Wigler, and J. Hicks. Optimizingsparse sequencing of single cells for highly multiplex copy numberprofiling. Genome Research, 125(5):714–724, 2015. → page 10[11] Y. Benjamini and T. P. Speed. Summarizing and correcting the GC contentbias in high-throughput sequencing. Nucleic acids research, 40(10):e72–e72, may 2012. → page 56[12] C. L. Bohrson, A. R. Barton, M. A. Lodato, R. E. Rodin, L. J. Luquette,V. V. Viswanadham, D. C. Gulhan, I. Cortés-Ciriano, M. A. Sherman,M. Kwon, M. E. Coulter, A. Galor, C. A. Walsh, and P. J. Park.Linked-read analysis identifies mutations in single-cell DNA-sequencingdata. Nature Genetics, 51(4):749–754, 2019. → pages 15, 58[13] D. D. Bowtell, S. Böhm, A. A. Ahmed, P.-J. Aspuria, R. C. Bast, V. Beral,J. S. Berek, M. J. Birrer, S. Blagden, M. A. Bookman, J. D. Brenton, K. B.Chiappinelli, F. C. Martins, G. Coukos, R. Drapkin, R. Edmondson,C. Fotopoulou, H. Gabra, J. Galon, C. Gourley, V. Heong, D. G. Huntsman,130M. Iwanicki, B. Y. Karlan, A. Kaye, E. Lengyel, D. A. Levine, K. H. Lu,I. A. McNeish, U. Menon, S. A. Narod, B. H. Nelson, K. P. Nephew,P. Pharoah, D. J. Powell, P. Ramos, I. L. Romero, C. L. Scott, A. K. Sood,E. A. Stronach, and F. R. Balkwill. Rethinking ovarian cancer II: reducingmortality from high-grade serous ovarian cancer. Nature Reviews Cancer,15(11):668–679, 2015. → page 82[14] I. Bozic, T. Antal, H. Ohtsuki, H. Carter, D. Kim, S. Chen, R. Karchin,K. W. Kinzler, B. Vogelstein, and M. A. Nowak. Accumulation of driverand passenger mutations during tumor progression. Proceedings of theNational Academy of Sciences, 107(43):18545–18550, 2010. → page 2[15] B. B. Campbell, N. Light, D. Fabrizio, M. Zatzman, F. Fuligni, R. de Borja,S. Davidson, M. Edwards, J. A. Elvin, K. P. Hodel, W. Zahurancik, Z. Suo,T. Lipman, K. Wimmer, C. P. Kratz, D. C. Bowers, T. W. Laetsch, G. P.Dunn, T. M. Johanns, M. R. Grimmer, I. V. Smirnov, V. Larouche,D. Samuel, A. Bronsema, M. Osborn, D. Stearns, P. Raman, K. A. Cole,P. B. Storm, M. Yalon, E. Opocher, G. Mason, G. A. Thomas, M. Sabel,B. George, D. S. Ziegler, S. Lindhorst, V. M. Issai, S. Constantini,H. Toledano, R. Elhasid, R. Farah, R. Dvir, P. Dirks, A. Huang, M. A.Galati, J. Chung, V. Ramaswamy, M. S. Irwin, M. Aronson, C. Durno,M. D. Taylor, G. Rechavi, J. M. Maris, E. Bouffet, C. Hawkins, J. F.Costello, M. S. Meyn, Z. F. Pursell, D. Malkin, U. Tabori, and A. Shlien.Comprehensive Analysis of Hypermutation in Human Cancer. Cell, 171:1042–1056.e10, 2017. → page 2[16] P. J. Campbell, S. Yachida, L. J. Mudie, P. J. Stephens, E. D. Pleasance,L. A. Stebbings, L. A. Morsberger, C. Latimer, S. McLaren, M.-L. Lin,D. J. McBride, I. Varela, S. A. Nik-Zainal, C. Leroy, M. Jia, A. Menzies,A. P. Butler, J. W. Teague, C. A. Griffin, J. Burton, H. Swerdlow, M. A.Quail, M. R. Stratton, C. Iacobuzio-Donahue, and P. A. Futreal. Thepatterns and dynamics of genomic instability in metastatic pancreaticcancer. Nature, 467(7319):1109–13, 2010. → page 3[17] B. L. Cantarel, D. Weaver, N. McNeill, J. Zhang, A. J. Mackey, andJ. Reese. Baysic: a bayesian method for combining sets of genome variantswith improved specificity and sensitivity. BMC bioinformatics, 15(1):104,2014. → page 14[18] S. Carreira, A. Romanel, J. Goodall, E. Grist, R. Ferraldeschi, S. Miranda,D. Prandi, D. Lorente, J.-S. Frenel, C. Pezaro, A. Omlin, D. N. Rodrigues,P. Flohr, N. Tunariu, J. S de Bono, F. Demichelis, and G. Attard. Tumor131clone dynamics in lethal prostate cancer. Science translational medicine, 6(254):254ra125, 2014. → page 3[19] J. Carrot-Zhang and J. Majewski. Lolopicker: detecting low allelic-fractionvariants from low-quality cancer samples. Oncotarget, 8(23):37032, 2017.→ pages 13, 14[20] Z. R. Chalmers, C. F. Connelly, D. Fabrizio, L. Gay, S. M. Ali, R. Ennis,A. Schrock, B. Campbell, A. Shlien, J. Chmielecki, F. Huang, Y. He,J. Sun, U. Tabori, M. Kennedy, D. S. Lieber, S. Roels, J. White, G. A. Otto,J. S. Ross, L. Garraway, V. A. Miller, P. J. Stephens, and G. M. Frampton.Analysis of 100,000 human cancer genomes reveals the landscape of tumormutational burden. Genome Medicine, 9(1):34, 2017. → page 2[21] D. Chen, H. Zhen, Y. Qiu, P. Liu, P. Zeng, J. Xia, Q. Shi, L. Xie, Z. Zhu,Y. Gao, G. Huang, J. Wang, H. Yang, and F. Chen. Comparison of singlecell sequencing data between two whole genome amplification methods ontwo sequencing platforms. Scientific Reports, 8(1):4963, 2018. → page 56[22] K. Cibulskis, M. S. Lawrence, S. L. Carter, A. Sivachenko, D. Jaffe,C. Sougnez, S. Gabriel, M. Meyerson, E. S. Lander, and G. Getz. Sensitivedetection of somatic point mutations in impure and heterogeneous cancersamples. Nature Biotechnology, 31(3):213–219, 2013. → pages 13, 20[23] P. Cingolani, A. Platts, L. L. Wang, M. Coon, T. Nguyen, L. Wang, S. J.Land, X. Lu, and D. M. Ruden. A program for annotating and predictingthe effects of single nucleotide polymorphisms, SnpEff: SNPs in thegenome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly, 6(2):80–92, 2012. → page 122[24] . G. P. Consortium. A map of human genome variation frompopulation-scale sequencing. Nature, 467(7319):1061–1073, 2010. →page 12[25] I. H. G. S. Consortium et al. Finishing the euchromatic sequence of thehuman genome. Nature, 431(7011):931, 2004. → page 6[26] C. S. Cooper, R. Eeles, D. C. Wedge, P. Van Loo, G. Gundem, L. B.Alexandrov, B. Kremeyer, A. Butler, A. G. Lynch, N. Camacho, C. E.Massie, J. Kay, H. J. Luxton, S. Edwards, Z. Kote-Jarai, N. Dennis,S. Merson, D. Leongamornlert, J. Zamora, C. Corbishley, S. Thomas,S. Nik-Zainal, S. O’Meara, L. Matthews, J. Clark, R. Hurst, R. Mithen,R. G. Bristow, P. C. Boutros, M. Fraser, S. Cooke, K. Raine, D. Jones,132A. Menzies, L. Stebbings, J. Hinton, J. Teague, S. McLaren, L. Mudie,C. Hardy, E. Anderson, O. Joseph, V. Goody, B. Robinson, M. Maddison,S. Gamble, C. Greenman, D. Berney, S. Hazell, N. Livni, C. Fisher,C. Ogden, P. Kumar, A. Thompson, C. Woodhouse, D. Nicol, E. Mayer,T. Dudderidge, N. C. Shah, V. Gnanapragasam, T. Voet, P. Campbell,A. Futreal, D. Easton, A. Y. Warren, C. S. Foster, M. R. Stratton, H. C.Whitaker, U. McDermott, D. S. Brewer, and D. E. Neal. Analysis of thegenetic phylogeny of multifocal prostate cancer identifies multipleindependent clonal expansions in neoplastic and morphologically normalprostate tissue. Nature genetics, 47(4):367–72, 2015. → page 3[27] E. C. de Bruin, N. McGranahan, R. Mitter, M. Salm, D. C. Wedge,L. Yates, M. Jamal-Hanjani, S. Shafi, N. Murugaesu, A. J. Rowan,E. Grönroos, M. A. Muhammad, S. Horswell, M. Gerlinger, I. Varela,D. Jones, J. Marshall, T. Voet, P. Van Loo, D. M. Rassl, R. C. Rintoul,S. M. Janes, S.-M. Lee, M. Forster, A. Tanya, D. Lawrence, M. Falzon,A. Capitanio, T. T. Harkins, C. C. Lee, W. Tom, E. Teefee, S.-C. Chen,S. Bengum, A. Rabinowitz, B. Phillimore, B. Spencer-Dene, G. Stamp,Z. Szallasi, N. Matthews, A. Stewart, P. Campbell, and C. Swanton. Spatialand temporal diversity in genomic instability processes defines lung cancerevolution. Science, 346(6206):251–256., 2014. → page 3[28] P. Del Moral, A. Doucet, and A. Jasra. Sequential monte carlo samplers.Journal of the Royal Statistical Society: Series B (Statistical Methodology),68(3):411–436, 2006. → page 64[29] M. G. Denis, A. Vallée, and S. Théoleyre. EGFR T790M resistancemutation in non small-cell lung carcinoma. Clinica Chimica Acta, 444:81–85, 2015. → page 3[30] A. G. Deshwar, S. Vembu, C. K. Yung, G. H. Jang, L. Stein, and Q. Morris.Phylowgs: reconstructing subclonal composition and evolution fromwhole-genome sequencing of tumors. Genome biology, 16(1):35, 2015. →page 19[31] J. Ding, A. Bashashati, A. Roth, A. Oloumi, K. Tse, T. Zeng, G. Haffari,M. Hirst, M. A. Marra, A. Condon, S. Aparicio, and S. P. Shah.Feature-based classifiers for somatic mutation detection in tumour-normalpaired sequencing data. Bioinformatics, 28(2):167–175, 2012. → pages14, 20, 35, 91133[32] L. Ding, M. J. Ellis, S. Li, D. E. Larson, K. Chen, J. W. Wallis, C. C.Harris, M. D. McLellan, R. S. Fulton, L. L. Fulton, et al. Genomeremodelling in a basal-like breast cancer metastasis and xenograft. Nature,464(7291):999, 2010. → pages 9, 10[33] X. Dong, L. Zhang, B. Milholland, M. Lee, A. Y. Maslov, T. Wang, andJ. Vijg. Accurate identification of single-nucleotide variants inwhole-genome-amplified single cells. Nature methods, 14(5):491—493,May 2017. → pages 15, 58[34] F. Dorri, S. Jewell, A. Bouchard-Côté, and S. P. Shah. Somatic mutationdetection and classification through probabilistic integration of clonalpopulation information. Communications Biology, 2(1):44, 2019. → page16[35] T. Dunn, G. Berry, D. Emig-Agius, Y. Jiang, S. Lei, A. Iyer, N. Udar, H.-Y.Chuang, J. Hegarty, M. Dickover, B. Klotzle, J. Robbins, M. Bibikova,M. Peeters, and M. Strömberg. Pisces: an accurate and versatile variantcaller for somatic and germline next-generation sequencing data.Bioinformatics, 35(9):1579–1581, 10 2018. → page 14[36] P. Eirew, A. Steif, J. Khattra, G. Ha, D. Yap, H. Farahani, K. Gelmon,S. Chia, C. Mar, A. Wan, et al. Dynamics of genomic clones in breastcancer patient xenografts at single-cell resolution. Nature, 518(7539):422,2015. → pages 3, 19[37] E. Falconer, M. Hills, U. Naumann, S. S. S. Poon, E. A. Chavez, A. D.Sanders, Y. Zhao, M. Hirst, and P. M. Lansdorp. DNA template strandsequencing of single-cells maps genomic rearrangements at highresolution. Nature methods, 9(11):1107–1112, nov 2012. → page 56[38] E. C. Friedberg, L. D. McDaniel, and R. A. Schultz. The role ofendogenous and exogenous DNA damage and mutagenesis. CurrentOpinion in Genetics & Development, 14(1):5–10, 2004. → page 2[39] E. Garrison and G. Marth. Haplotype-based variant detection fromshort-read sequencing. Preprint at arXiv:1207.3907v2, 2012. → pages13, 20, 35[40] C. Gawad, W. Koh, and S. R. Quake. Dissecting the clonal origins ofchildhood acute lymphoblastic leukemia by single-cell genomics.Proceedings of the National Academy of Sciences of the United States ofAmerica, 111(50):17947–52, 2014. → page 56134[41] C. Gawad, W. Koh, and S. R. Quake. Single-cell genome sequencing:current state of the science. Nature Review Genetics, 17(3):175–188, 2016.→ page 56[42] M. Gerlinger, A. J. Rowan, S. Horswell, J. Larkin, D. Endesfelder,E. Gronroos, P. Martinez, N. Matthews, A. Stewart, P. Tarpey, I. Varela,B. Phillimore, S. Begum, N. Q. McDonald, A. Butler, D. Jones, K. Raine,C. Latimer, C. R. Santos, M. Nohadani, A. C. Eklund, B. Spencer-Dene,G. Clark, L. Pickering, G. Stamp, M. Gore, Z. Szallasi, J. Downward, P. A.Futreal, and C. Swanton. Intratumor Heterogeneity and BranchedEvolution Revealed by Multiregion Sequencing. 366(10):883–892, 2012.→ pages 3, 19[43] M. Gerlinger, S. Horswell, J. Larkin, A. J. Rowan, M. P. Salm, I. Varela,R. Fisher, N. McGranahan, N. Matthews, C. R. Santos, P. Martinez,B. Phillimore, S. Begum, A. Rabinowitz, B. Spencer-Dene, S. Gulati, P. A.Bates, G. Stamp, L. Pickering, M. Gore, D. L. Nicol, S. Hazell, P. A.Futreal, A. Stewart, and C. Swanton. Genomic architecture and evolutionof clear cell renal cell carcinomas defined by multiregion sequencing.Nature genetics, 46(3):225–33, 2014. → page 3[44] M. Gerstung, C. Beisel, M. Rechsteiner, P. Wild, P. Schraml, H. Moch, andN. Beerenwinkel. Reliable detection of subclonal single-nucleotidevariants in tumour cell populations. Nature communications, 3:811, 2012.→ pages 13, 14[45] O. Gevaert, V. Villalobos, B. I. Sikic, and S. K. Plevritis. Identification ofovarian cancer driver genes by using module network integration ofmulti-omics data. Interface focus, 3(4):20130013, aug 2013. → page 82[46] R. Govindan, L. Ding, M. Griffith, J. Subramanian, N. D. Dees, K. L.Kanchi, C. A. Maher, R. Fulton, L. Fulton, J. Wallis, K. Chen, J. Walker,S. Mcdonald, R. Bos, D. Ornitz, D. Xiong, M. You, D. J. Dooling, andM. Watson. GENOMIC LANDSCAPE OF NON-SMALL CELL LUNGCANCER IN SMOKERS AND NEVER SMOKERS. Cell, 150(6):1121–1134, 2012. → page 19[47] R. Goya, M. G. F. Sun, R. D. Morin, G. Leung, G. Ha, K. C. Wiegand,J. Senz, A. Crisan, M. A. Marra, M. Hirst, D. Huntsman, K. P. Murphy,S. Aparicio, and S. P. Shah. SNVMix: Predicting single nucleotide variantsfrom next-generation sequencing of tumors. Bioinformatics, 26(6):730–736, 2010. → pages 14, 20135[48] M. Greaves and C. C. Maley. Clonal evolution in cancer. Nature, 481(7381):306, 2012. → page 2[49] W. W. R. H. J. G. Bloom. Histological grading and prognosis of breastcancer. British Journal of Cancer, 11(3):359–377, 1957. → page 3[50] G. Ha, A. Roth, J. Khattra, J. Ho, D. Yap, L. M. Prentice, N. Melnyk,A. McPherson, A. Bashashati, E. Laks, et al. Titan: inference of copynumber architectures in clonal cell populations from tumor whole-genomesequence data. Genome research, 24(11):1881–1893, 2014. → pages21, 35[51] M. C. Haffner, T. Mosbruger, D. M. Esopi, H. Fedor, C. M. Heaphy, D. A.Walker, N. Adejola, M. Gurel, J. Hicks, A. K. Meeker, M. K. Halushka,J. W. Simons, W. B. Isaacs, A. M. De Marzo, W. G. Nelson, andS. Yegnasubramanian. Tracking the clonal origin of lethal prostate cancer.Journal of Clinical Investigation, 123(11):4918–4922, 2013. → page 3[52] D. Hanahan and R. A. Weinberg. The hallmarks of cancer. cell, 100(1):57–70, 2000. → page 2[53] N. F. Hansen, J. J. Gartner, L. Mei, Y. Samuels, and J. C. Mullikin.Shimmer: detection of genetic alterations in tumors using next-generationsequence data. Bioinformatics, 29(12):1498–1503, 2013. → page 14[54] J. Hård, E. Al Hakim, M. Kindblom, Å. K. Björklund, B. Sennblad,I. Demirci, M. Paterlini, P. Reu, E. Borgström, P. L. Ståhl, J. Michaelsson,J. E. Mold, and J. Frisén. Conbase: a software for unsupervised discoveryof clonal somatic mutations in single cells through read phasing. 20(1):68,2019. → pages 15, 58[55] G. H. Heppner and B. E. Miller. Tumor heterogeneity: biologicalimplications and therapeutic consequences. Cancer and Metastasis Review,2(1):5–23, 1983. → page 3[56] Y. Hou, L. Song, P. Zhu, B. Zhang, Y. Tao, X. Xu, F. Li, K. Wu, J. Liang,D. Shao, H. Wu, X. Ye, C. Ye, R. Wu, M. Jian, Y. Chen, W. Xie, R. Zhang,L. Chen, X. Liu, X. Yao, H. Zheng, C. Yu, Q. Li, Z. Gong, M. Mao,X. Yang, L. Yang, J. Li, W. Wang, Z. Lu, N. Gu, G. Laurie, L. Bolund,K. Kristiansen, J. Wang, H. Yang, Y. Li, X. Zhang, and J. Wang.Single-Cell Exome Sequencing and Monoclonal Evolution of aJAK2-Negative Myeloproliferative Neoplasm. Cell, 148(5):873–885, 2012.→ page 56136[57] X. Huang, X.-n. Wang, X.-d. Yuan, W.-y. Wu, P. E. Lobie, and Z. Wu.XIAP facilitates breast and colon carcinoma growth via promotion of p62depletion through ubiquitination-dependent proteasomal degradation.Oncogene, 38(9):1448–1460, 2019. → page 122[58] S. C. Huhn, J. Liu, C. Ye, H. Lu, X. Jiang, X. Feng, S. Ganesan, E. White,and Z. Shen. Regulation of spindle integrity and mitotic fidelity by BCCIP.Oncogene, 36(33):4750–4766, 2017. → page 122[59] Illumina. An introduction to next generation sequencing technology,January 2011. URL http://www.illumina.com/content/dam/illumina-marketing/documents/products/illumina_sequencing_introduction.pdf. →page 8[60] M. Jamal-Hanjani, G. A. Wilson, N. McGranahan, N. J. Birkbak, T. B.Watkins, S. Veeriah, S. Shafi, D. H. Johnson, R. Mitter, R. Rosenthal, et al.Tracking the evolution of non–small-cell lung cancer. New EnglandJournal of Medicine, 376(22):2109–2121, 2017. → pages 19, 37, 38, 40[61] R. Jeselsohn, G. Buchwalter, C. De Angelis, M. Brown, and R. Schiff.ESR1 mutations—a mechanism for acquired endocrine resistance in breastcancer. Nature reviews. Clinical oncology, 12(10):573–583, oct 2015. →page 3[62] Q. Jin, X. Mao, B. Li, S. Guan, F. Yao, and F. Jin. Overexpression ofSMARCA5 correlates with cell proliferation and migration in breastcancer. Tumor Biology, 36(3):1895–1902, 2015. → page 122[63] M. Josephidou, A. G. Lynch, and S. Tavaré. multiSNV: a probabilisticapproach for improving detection of somatic point mutations from multiplerelated tumour samples. Nucleic acids research, 43(9):e61, 2015. → pages15, 20, 35[64] I. Kalatskaya, Q. M. Trinh, M. Spears, J. D. McPherson, J. M. Bartlett, andL. Stein. Isown: accurate somatic mutation identification in the absence ofnormal tissue controls. Genome medicine, 9(1):59, 2017. → pages 14, 15[65] Z. Kan, B. S. Jaiswal, J. Stinson, V. Janakiraman, D. Bhatt, H. M. Stern,P. Yue, P. M. Haverty, R. Bourgon, J. Zheng, M. Moorhead, S. Chaudhuri,L. P. Tomsho, B. A. Peters, K. Pujara, S. Cordes, D. P. Davis, V. E. H.Carlton, W. Yuan, L. Li, W. Wang, C. Eigenbrot, J. S. Kaminker, D. A.Eberhard, P. Waring, S. C. Schuster, Z. Modrusan, Z. Zhang, D. Stokoe,F. J. de Sauvage, M. Faham, and S. Seshagiri. Diverse somatic mutation137patterns and pathway alterations in human cancers. Nature, 466(7308):869–873, 2010. → page 82[66] C. Kandoth, M. D. McLellan, F. Vandin, K. Ye, B. Niu, C. Lu, M. Xie,Q. Zhang, J. F. McMichael, M. A. Wyczalkowski, M. D. M. Leiserson,C. A. Miller, J. S. Welch, M. J. Walter, M. C. Wendl, T. J. Ley, R. K.Wilson, B. J. Raphael, and L. Ding. Mutational landscape and significanceacross 12 major cancer types. Nature, 502(7471):333–339, 2013. → page 3[67] S. Kim, K. Jeong, K. Bhutani, J. H. Lee, A. Patel, E. Scott, H. Nam,H. Lee, J. G. Gleeson, and V. Bafna. Virmid: accurate detection of somaticmutations with sample impurity inference. Genome biology, 14(8):R90,2013. → pages 12, 13, 20[68] M. Kleppe and R. L. Levine. Tumor heterogeneity confounds andilluminates: assessing the implications. Nature medicine, 20(4):342–4,2014. → page 3[69] D. C. Koboldt, K. Chen, T. Wylie, D. E. Larson, M. D. McLellan, E. R.Mardis, G. M. Weinstock, R. K. Wilson, and L. Ding. VarScan: Variantdetection in massively parallel sequencing of individual and pooledsamples. Bioinformatics, 25(17):2283–2285, 2009. → pages 14, 20[70] C. Kockan, F. Hach, I. Sarrafi, R. H. Bell, B. McConeghy, K. Beja,A. Haegert, A. W. Wyatt, S. V. Volik, K. N. Chi, et al. Sinvict:ultra-sensitive detection of single nucleotide variants and indels incirculating tumour dna. Bioinformatics, 33(1):26–34, 2016. → pages14, 15[71] R. Kridel, F. C. Chan, A. Mottok, M. Boyle, P. Farinha, K. Tan,B. Meissner, A. Bashashati, A. McPherson, A. Roth, et al. Histologicaltransformation and progression in follicular lymphoma: a clonal evolutionstudy. PLoS medicine, 13(12):e1002197, 2016. → page 19[72] J. Kuipers, P. Suter, and G. Moffa. Efficient structure learning andsampling of bayesian networks. 2018. → pages 57, 77[73] S. Kumar and S. Subramanian. Mutation rates in mammalian genomes.Proceedings of the National Academy of Sciences, 99(2):803–808, 2002.→ page 2[74] Z. Lai, A. Markovets, M. Ahdesmaki, B. Chapman, O. Hofmann,R. McEwen, J. Johnson, B. Dougherty, J. C. Barrett, and J. R. Dry. Vardict:138a novel and versatile variant caller for next-generation sequencing in cancerresearch. Nucleic acids research, 44(11):e108–e108, 2016. → page 14[75] E. Laks, H. Zahn, D. Lai, A. McPherson, A. Steif, J. Brimhall, J. Biele,B. Wang, T. Masud, D. Grewal, C. Nielsen, S. Leung, V. Bojilova,M. Smith, O. Golovko, S. Poon, P. Eirew, F. Kabeer, T. R. de Algara, S. R.Lee, M. J. Taghiyar, C. Huebner, J. Ngo, T. Chan, S. Vatrt-Watts,P. Walters, N. Abrar, S. Chan, M. Wiens, L. Martin, R. W. Scott, M. T.Underhill, E. Chavez, C. Steidl, D. D. Costa, Y. Ma, R. J. N. Coope,R. Corbett, S. Pleasance, R. Moore, A. J. Mungall, M. A. Marra,C. Hansen, S. P. Shah, and S. Aparicio. Resource: Scalable whole genomesequencing of 40,000 single cells identifies stochastic aneuploidies,genome replication states and clonal repertoires. bioRxiv, 2018. → page 56[76] E. Laks, A. McPherson, H. Zahn, D. Lai, A. Steif, J. Brimhall, J. Biele,B. Wang, T. Masud, J. Ting, D. Grewal, C. Nielsen, S. Leung, V. Bojilova,M. Smith, O. Golovko, S. Poon, P. Eirew, F. Kabeer, T. Ruiz de Algara,S. R. Lee, M. J. Taghiyar, C. Huebner, J. Ngo, T. Chan, S. Vatrt-Watts,P. Walters, N. Abrar, S. Chan, M. Wiens, L. Martin, R. W. Scott, T. M.Underhill, E. Chavez, C. Steidl, D. Da Costa, Y. Ma, R. J. N. Coope,R. Corbett, S. Pleasance, R. Moore, A. J. Mungall, C. Mar, F. Cafferty,K. Gelmon, S. Chia, G. J. Hannon, G. Battistoni, D. Bressan, I. Cannell,H. Casbolt, C. Jauset, T. Kovacˇevic´, C. Mulvey, F. Nugent, M. P. Ribes,I. Pearsall, F. Qosaj, K. Sawicka, S. Wild, E. Williams, S. Aparicio,E. Laks, Y. Li, C. O’Flanagan, A. Smith, T. Ruiz, S. Balasubramanian,M. Lee, B. Bodenmiller, M. Burger, L. Kuett, S. Tietscher, J. Windager,E. Boyden, S. Alon, Y. Cui, A. Emenari, D. Goodwin, E. Karagiannis,A. Sinha, A. T. Wassie, C. Caldas, A. Bruna, M. Callari, W. Greenwood,G. Lerda, Y. Lubling, A. Marti, O. Rueda, A. Shea, O. Harris, R. Becker,F. Grimaldi, S. Harris, S. Vogl, J. A. Joyce, J. Hausser, S. Watson, S. Shah,A. McPherson, I. Vázquez-García, S. Tavaré, K. Dinh, E. Fisher, R. Kunes,N. A. Walton, M. Al Sa’d, N. Chornay, A. Dariush, E. G. Solares,C. Gonzalez-Fernandez, A. K. Yoldas, N. Millar, X. Zhuang, J. Fan, H. Lee,L. S. Duran, C. Xia, P. Zheng, M. A. Marra, C. Hansen, S. P. Shah, andS. Aparicio. Clonal Decomposition and DNA Replication States Definedby Scaled Single-Cell Genome Sequencing. Cell, 179(5):1207–1221.e22,2019. → pages 11, 82, 84, 88, 95, 96, 102, 103, 104, 105, 116, 117[77] D. A. Landau, S. L. Carter, P. Stojanov, A. McKenna, K. Stevenson, M. S.Lawrence, C. Sougnez, C. Stewart, A. Sivachenko, L. Wang, Y. Wan,W. Zhang, S. A. Shukla, A. Vartanov, S. M. Fernandes, G. Saksena,139K. Cibulskis, B. Tesar, S. Gabriel, N. Hacohen, M. Meyerson, E. S. Lander,D. Neuberg, J. R. Brown, G. Getz, and C. J. Wu. Evolution and impact ofsubclonal mutations in chronic lymphocytic leukemia. Cell, 152(4):714–726, 2013. → page 3[78] E. S. Lander, L. M. Linton, B. Birren, C. Nusbaum, M. C. Zody,J. Baldwin, K. Devon, K. Dewar, M. Doyle, W. Fitzhugh, et al. Erratum:Initial sequencing and analysis of the human genome: International humangenome sequencing consortium (nature (2001) 409 (860-921)). Nature,412(6846):565–566, 2001. → page 6[79] D. Larson, C. Harris, K. Chen, D. Koboldt, T. Abbott, D. Dooling, T. Ley,E. Mardis, R. Wilson, and L. Ding. Somaticsniper: identification ofsomatic point mutations in whole genome sequencing data. Bioinformatics(Oxford, England), 28:311–7, 12 2011. → page 12[80] A. Lex, N. Gehlenborg, H. Strobelt, R. Vuillemot, and H. Pfister. Upset:visualization of intersecting sets. IEEE transactions on visualization andcomputer graphics, 20(12):1983–1992, 2014. → page 53[81] Y. Liu, M. Loewer, S. Aluru, and B. Schmidt. Snvsniffer: an integratedcaller for germline and somatic single-nucleotide and indel mutations.BMC systems biology, 10(2):47, 2016. → pages 12, 13[82] I. C. Macaulay and T. Voet. Single Cell Genomics: Advances and FuturePerspectives. PLoS Genetics, 10(1):1–9, 2014. → page 11[83] A. McKenna, M. Hanna, E. Banks, A. Sivachenko, K. Cibulskis,A. Kernytsky, K. Garimella, D. Altshuler, S. Gabriel, M. Daly, et al. Thegenome analysis toolkit: a mapreduce framework for analyzingnext-generation dna sequencing data. Genome research, 20(9):1297–1303,2010. → pages 15, 20[84] A. McPherson, A. Roth, E. Laks, T. Masud, A. Bashashati, A. W. Zhang,G. Ha, J. Biele, D. Yap, A. Wan, L. M. Prentice, J. Khattra, M. A. Smith,C. B. Nielsen, S. C. Mullaly, S. Kalloger, A. Karnezis, K. Shumansky,C. Siu, J. Rosner, H. L. Chan, J. Ho, N. Melnyk, J. Senz, W. Yang,R. Moore, A. J. Mungall, M. A. Marra, A. Bouchard-Côté, C. B. Gilks,D. G. Huntsman, J. N. McAlpine, S. Aparicio, and S. P. Shah. Divergentmodes of clonal spread and intraperitoneal mixing in high-grade serousovarian cancer. Nature Genetics, 48(October 2015):758–767, 2016. →pages 19, 34, 36, 37140[85] X. Meng, J. Fan, and Z. Shen. Roles of BCCIP in chromosome stabilityand cytokinesis. Oncogene, 26(43):6253–6260, sep 2007. → page 122[86] E. Muller, N. Goardon, B. Brault, A. Rousselin, G. Paimparay, A. Legros,R. Fouillet, O. Bruet, A. Tranchant, F. Domin, et al. Outlyzer: software forextracting low-allele-frequency tumor mutations from sequencingbackground noise in clinical practice. Oncotarget, 7(48):79485, 2016. →page 14[87] G. H. Murillo, N. You, X. Su, W. Cui, M. P. Reilly, M. Li, K. Ning, andX. Cui. MultiGeMS: detection of SNVs from multiple samples usingmodel selection on high-throughput sequencing data. Bioinformatics, 32(10):1486–1492, 01 2016. → page 15[88] N. Navin, A. Krasnitz, L. Rodgers, K. Cook, J. Meth, J. Kendall, M. Riggs,Y. Eberling, J. Troge, V. Grubor, D. Levy, P. Lundin, S. Maner,A. Zetterberg, J. Hicks, and M. Wigler. Inferring tumor progression fromgenomic heterogeneity. Genome Research, 20(1):68–80, 2010. → page 2[89] N. Navin, J. Kendall, J. Troge, P. Andrews, L. Rodgers, J. McIndoo,K. Cook, A. Stepansky, D. Levy, D. Esposito, L. Muthuswamy,A. Krasnitz, W. R. McCombie, J. Hicks, and M. Wigler. Tumour evolutioninferred by single-cell sequencing. Nature, 472(7341):90–94, 2011. →pages 3, 10, 56, 60[90] S. Nik-Zainal, P. Van Loo, D. C. Wedge, L. B. Alexandrov, C. D.Greenman, K. W. Lau, K. Raine, D. Jones, J. Marshall, M. Ramakrishna,A. Shlien, S. L. Cooke, J. Hinton, A. Menzies, L. A. Stebbings, C. Leroy,M. Jia, R. Rance, L. J. Mudie, S. J. Gamble, P. J. Stephens, S. McLaren,P. S. Tarpey, E. Papaemmanuil, H. R. Davies, I. Varela, D. J. McBride,G. R. Bignell, K. Leung, A. P. Butler, J. W. Teague, S. Martin, G. Jonsson,O. Mariani, S. Boyault, P. Miron, A. Fatima, A. Langerod, S. A. J. R.Aparicio, A. Tutt, A. M. Sieuwerts, k. Borg, G. Thomas, A. V. Salomon,A. L. Richardson, A. L. Borresen-Dale, P. A. Futreal, M. R. Stratton, andP. J. Campbell. The life history of 21 breast cancers. Cell, 149(5):994–1007, 2012. → pages 3, 19[91] S. M. Noordermeer and H. van Attikum. PARP Inhibitor Resistance: ATug-of-War in BRCA-Mutated Cells. Trends in Cell Biology, 29(10):820–834, oct 2019. → page 3[92] P. C. Nowell. The clonal evolution of tumor cell populations. Science, 194(4260):23–28, 1988. → page 2141[93] S. E. Park, K. Park, E. Lee, J.-Y. Kim, J. S. Ahn, Y.-H. Im, C. Lee, H. Jung,S. Y. Cho, W.-Y. Park, R. Cristescu, and Y. H. Park. Clinical implication oftumor mutational burden in patients with HER2-positive refractorymetastatic breast cancer. OncoImmunology, 7(8):e1466768, aug 2018. →page 82[94] E. D. Pleasance, P. J. Stephens, S. O’Meara, D. J. McBride, A. Meynert,D. Jones, M.-L. Lin, D. Beare, K. W. Lau, C. Greenman, et al. A small-celllung cancer genome with complex signatures of tobacco exposure. Nature,463(7278):184, 2010. → page 9[95] V. Popic, R. Salari, I. Hajirasouliha, D. Kashef-Haghighi, R. B. West, andS. Batzoglou. Fast and scalable inference of multi-sample cancer lineages.Genome biology, 16(1):91, 2015. → page 19[96] P. Priestley, J. Baber, M. P. Lolkema, N. Steeghs, E. de Bruijn, C. Shale,K. Duyvesteyn, S. Haidari, A. van Hoeck, W. Onstenk, P. Roepman,M. Voda, H. J. Bloemendal, V. C. Tjan-Heijnen, C. M. van Herpen,M. Labots, P. O. Witteveen, E. F. Smit, S. Sleijfer, E. E. Voest, andE. Cuppen. Pan-cancer whole genome analyses of metastatic solid tumors.bioRxiv, 2019. → page 9[97] B. K. Rajendran and C.-X. Deng. Characterization of potential drivermutations involved in human breast cancer by computational approaches.Oncotarget, 8(30):50252–50272, July 2017. → page 82[98] A. Rimmer, H. Phan, I. Mathieson, Z. Iqbal, S. R. F. Twigg, A. O. M.Wilkie, G. McVean, and G. Lunter. Integrating mapping-, assembly- andhaplotype-based approaches for calling variants in clinical sequencingapplications. Nature genetics, 46(8):912–918, 2014. → page 13[99] A. Roth, J. Ding, R. Morin, A. Crisan, G. Ha, R. Giuliany, A. Bashashati,M. Hirst, G. Turashvili, A. Oloumi, M. A. Marra, S. Aparicio, and S. P.Shah. JointSNVMix: A probabilistic model for accurate detection ofsomatic mutations in normal/tumour paired next-generation sequencingdata. Bioinformatics, 28(7):907–913, 2012. → pages 12, 13, 20[100] A. Roth, J. Khattra, D. Yap, A. Wan, E. Laks, J. Biele, G. Ha, S. Aparicio,A. Bouchard-Côté, and S. P. Shah. PyClone: statistical inference of clonalpopulation structure in cancer. Nature Methods, 11(4):396–398, 2014. →pages 9, 19, 34142[101] A. Roth, A. McPherson, E. Laks, J. Biele, D. Yap, A. Wan, M. A. Smith,C. B. Nielsen, J. N. McAlpine, S. Aparicio, A. Bouchard-Côté, and S. P.Shah. Clonal genotype and population structure inference from single-celltumor sequencing. Nature Methods, 2016. → page 16[102] A. J. L. Roth. Probabilistic models for the identification and interpretationof somatic single nucleotide variants in cancer genomes. PhD thesis,University of British Columbia, 2015. → page 56[103] R. Salari, S. S. Saleh, D. Kashef-Haghighi, D. Khavari, D. E. Newburger,R. B. West, A. Sidow, and S. Batzoglou. Inference of tumor phylogenieswith improved somatic mutation discovery. Journal of ComputationalBiology, 20(11):933–944, 2013. → page 20[104] S. Salehi, F. Kabir, A. Bouchard-Côté, S. P. Shah, et al. Quantitativegenotype and phenotype fitness modeling in cancerusing single celltimeseries population dynamics. → pages84, 89, 90, 98, 99, 100, 101, 107, 108, 109, 110, 111, 112, 113, 114, 118, 119, 120, 121[105] C. T. Saunders, W. S. W. Wong, S. Swamy, J. Becq, L. J. Murray, and R. K.Cheetham. Strelka: Accurate somatic small-variant calling from sequencedtumor-normal sample pairs. Bioinformatics, 28(14):1811–1817, 2012. →pages 13, 20, 35, 39, 91[106] R. F. Schwarz, C. K. Y. Ng, S. L. Cooke, S. Newman, J. Temple, A. M.Piskorz, D. Gale, K. Sayal, M. Murtaza, P. J. Baldwin, N. Rosenfeld, H. M.Earl, E. Sala, M. Jimenez-Linan, C. A. Parkinson, F. Markowetz, and J. D.Brenton. Spatial and Temporal Heterogeneity in High-Grade SerousOvarian Cancer: A Phylogenetic Analysis. PLoS Medicine, 12(2):1–20,2015. → page 3[107] S. P. Shah, R. D. Morin, J. Khattra, L. Prentice, T. Pugh, A. Burleigh,A. Delaney, K. Gelmon, R. Guliany, J. Senz, C. Steidl, R. A. Holt, S. Jones,M. Sun, G. Leung, R. Moore, T. Severson, G. A. Taylor, A. E.Teschendorff, K. Tse, G. Turashvili, R. Varhol, R. L. Warren, P. Watson,Y. Zhao, C. Caldas, D. Huntsman, M. Hirst, M. A. Marra, and S. Aparicio.Mutational evolution in a lobular breast tumour profiled at single nucleotideresolution. Nature, 461(7265):809–813, 2009. → pages 9, 10, 19[108] 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,143R. 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. G. 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. A. Marra, and S. Aparicio. Theclonal and mutational evolution spectrum of primary triple-negative breastcancers. Nature, 486(7403):395–9, 2012. → pages 3, 9, 10[109] E. Shapiro, T. Biezuner, and S. Linnarsson. Single-cell sequencing-basedtechnologies will revolutionize whole-organism science. Nature ReviewsGenetics, 14:618, jul 2013. → page 56[110] G. N. Sharma, R. Dave, J. Sanadya, P. Sharma, and K. K. Sharma. Varioustypes and management of breast cancer: an overview. Journal of advancedpharmaceutical technology & research, 1(2):109–126, apr 2010. → page 83[111] Y. Shiraishi, Y. Sato, K. Chiba, Y. Okuno, Y. Nagata, K. Yoshida, N. Shiba,Y. Hayashi, H. Kume, Y. Homma, M. Sanada, S. Ogawa, and S. Miyano.An empirical Bayesian framework for somatic mutation detection fromcancer genome sequencing data. Nucleic Acids Research, 41(7):e89–e89,03 2013. → page 13[112] A. D. Siddiqui and B. Piperdi. KRAS Mutation in Colon Cancer: A Markerof Resistance to EGFR-I Therapy. Annals of Surgical Oncology, 17(4):1168–1176, 2010. → page 3[113] M. B. Siegel, X. He, K. A. Hoadley, A. Hoyle, J. B. Pearce, A. L. Garrett,S. Kumar, V. J. Moylan, C. M. Brady, A. E. Van Swearingen, et al.Integrated rna and dna sequencing reveals early drivers of metastatic breastcancer. The Journal of clinical investigation, 128(4), 2018. → page 19[114] N. Silver, S. Best, J. Jiang, and S. L. Thein. Selection of housekeepinggenes for gene expression studies in human reticulocytes using real-timePCR. BMC molecular biology, 7:33, oct 2006. → page 82[115] J. Singer, J. Kuipers, K. Jahn, and N. Beerenwinkel. Single-cell mutationidentification via phylogenetic inference. Nature Communications, 9(1):5144, 2018. → pages 16, 59144[116] K. S. Smith, V. K. Yadav, S. Pei, D. A. Pollyea, C. T. Jordan, and S. De.Somvarius: somatic variant identification from unpaired tissue samples.Bioinformatics, 32(6):808–813, 2015. → pages 14, 15[117] A. Sottoriva, I. Spiteri, S. G. M. Piccirillo, A. Touloumis, V. P. Collins,J. C. Marioni, C. Curtis, C. Watts, and S. Tavaré. Intratumor heterogeneityin human glioblastoma reflects cancer evolutionary dynamics. Proceedingsof the National Academy of Sciences of the United States of America, 110(10):4009–14, 2013. → page 3[118] A. Sottoriva, H. Kang, Z. Ma, T. A. Graham, M. P. Salomon, J. Zhao,P. Marjoram, K. Siegmund, M. F. Press, D. Shibata, and C. Curtis. A BigBang model of human colorectal tumor growth. Nature Genetics, 47(3):209–216, 2015. → page 3[119] J.-F. Spinella, P. Mehanna, R. Vidal, V. Saillour, P. Cassart, C. Richer,M. Ouimet, J. Healy, and D. Sinnett. Snooper: a machine learning-basedmethod for somatic variant identification from low-pass next-generationsequencing. BMC genomics, 17(1):912, 2016. → page 14[120] M. R. Stratton. Exploring the genomes of cancer cells: progress andpromise. science, 331(6024):1553–1558, 2011. → page 2[121] H. Suzuki, K. Aoki, K. Chiba, Y. Sato, Y. Shiozawa, Y. Shiraishi,T. Shimamura, A. Niida, K. Motomura, F. Ohka, T. Yamamoto,K. Tanahashi, M. Ranjit, T. Wakabayashi, T. Yoshizato, K. Kataoka,K. Yoshida, Y. Nagata, A. Sato-Otsubo, H. Tanaka, M. Sanada, Y. Kondo,H. Nakamura, M. Mizoguchi, T. Abe, Y. Muragaki, R. Watanabe, I. Ito,S. Miyano, A. Natsume, and S. Ogawa. Mutational landscape and clonalarchitecture in grade II and III gliomas. Nat Genet, 47(5):458–468, 2015.→ page 3[122] C. Thirlwell, O. C. C. Will, E. Domingo, T. A. Graham, S. A. C.McDonald, D. Oukrif, R. Jeffrey, M. Gorman, M. Rodriguez-Justo,J. Chin-Aleong, S. K. Clark, M. R. Novelli, J. A. Jankowski, N. A. Wright,I. P. M. Tomlinson, and S. J. Leedham. Clonality Assessment and ClonalOrdering of Individual Neoplastic Crypts Shows Polyclonality ofColorectal Adenomas. Gastroenterology, 12(11):1058–1060, 2010. →page 3[123] A. a. Tone, M. K. McConechy, W. Yang, J. Ding, S. Yip, E. Kong, K.-K.Wong, D. M. Gershenson, H. Mackay, S. Shah, B. Gilks, A. V. Tinker,145B. Clarke, J. N. McAlpine, and D. Huntsman. Intratumoral heterogeneityin a minority of ovarian low-grade serous carcinomas. BMC cancer, 14(1):982, 2014. → page 3[124] F. Vallania, E. Ramos, S. Cresci, R. D. Mitra, and T. E. Druley. Detectionof rare genomic variants from pooled sequencing using splinter. JoVE(Journal of Visualized Experiments), (64):e3943, 2012. → page 14[125] P. Van Loo, S. H. Nordgard, O. C. Lingjærde, H. G. Russnes, I. H. Rye,W. Sun, V. J. Weigman, P. Marynen, A. Zetterberg, B. Naume, C. M. Perou,A.-L. Børresen-Dale, and V. N. Kristensen. Allele-specific copy numberanalysis of tumors. Proceedings of the National Academy of Sciences, 107(39):16910–16915, 2010. → page 21[126] K. E. Van Rens, V. Mäkinen, and A. I. Tomescu. SNV-PPILP: RefinedSNV calling for tumor data using perfect phylogenies and ILP.Bioinformatics, 31(7):1133–1135, 2015. → pages 15, 20[127] W. Wang, P. Wang, F. Xu, R. Luo, M. P. Wong, T.-W. Lam, and J. Wang.Fasd-somatic: a fast and accurate somatic snv detection algorithm forcancer genome sequencing data. Bioinformatics, 30(17):2498–2500, 2014.→ page 12[128] Y. Wang, J. Waters, M. L. Leung, A. Unruh, W. Roh, X. Shi, K. Chen,P. Scheet, S. Vattathil, H. Liang, A. Multani, H. Zhang, R. Zhao, F. Michor,F. Meric-Bernstam, and N. E. Navin. Clonal evolution in breast cancerrevealed by single nucleus genome sequencing. Nature, 512(7513):1–15,2014. → pages 10, 56, 60[129] Z. Wei, W. Wang, P. Hu, G. J. Lyon, and H. Hakonarson. Snver: astatistical tool for variant calling in analysis of pooled or individualnext-generation sequencing data. Nucleic acids research, 39(19):e132–e132, 2011. → page 14[130] M. J. Williams, A. Sottoriva, and T. A. Graham. Measuring clonalevolution in cancer with genomics. Annual Review of Genomics andHuman Genetics, 20(1):309–329, 2019. → page 2[131] A. Wilm, P. P. K. Aw, D. Bertrand, G. H. T. Yeo, S. H. Ong, C. H. Wong,C. C. Khor, R. Petric, M. L. Hibberd, and N. Nagarajan. Lofreq: asequence-quality aware, ultra-sensitive variant caller for uncoveringcell-population heterogeneity from high-throughput sequencing datasets.Nucleic acids research, 40(22):11189–11201, 2012. → page 13146[132] L. D. Wood, D. W. Parsons, S. Jones, J. Lin, T. Sjöblom, R. J. Leary,D. Shen, S. M. Boca, T. Barber, J. Ptak, N. Silliman, S. Szabo, Z. Dezso,V. Ustyanksky, T. Nikolskaya, Y. Nikolsky, R. Karchin, P. A. Wilson, J. S.Kaminker, Z. Zhang, R. Croshaw, J. Willis, D. Dawson, M. Shipitsin,J. K. V. Willson, S. Sukumar, K. Polyak, B. H. Park, C. L. Pethiyagoda,P. V. K. Pant, D. G. Ballinger, A. B. Sparks, J. Hartigan, D. R. Smith,E. Suh, N. Papadopoulos, P. Buckhaults, S. D. Markowitz, G. Parmigiani,K. W. Kinzler, V. E. Velculescu, and B. Vogelstein. The genomiclandscapes of human breast and colorectal cancers. Science, 318(5853):1108–1113, 2007. → page 82[133] F. Wu. Irreversible Parallel Tempering and an Application to a BayesianNonparametric. PhD thesis, Oxford University, 2017. → pages 64, 73[134] T. Xie, M. Musteanu, P. P. Lopez-Casas, D. J. Shields, P. Olson, P. A. Rejto,and M. Hidalgo. Whole exome sequencing of rapid autopsy tumors andxenograft models reveals possible driver mutations underlying tumorprogression. PloS one, 10(11):e0142631, 2015. → page 19[135] C. Xu. A review of somatic single nucleotide variant calling algorithms fornext-generation sequencing data. Computational and StructuralBiotechnology Journal, 16:15–24, 2018. → pages 11, 12, 13, 14[136] S. Yachida, S. Jones, I. Bozic, T. Antal, R. Leary, B. Fu, M. Kamiyama,R. H. Hruban, J. R. Eshleman, M. A. Nowak, V. E. Velculescu, K. W.Kinzler, B. Vogelstein, and C. A. Iacobuzio-Donahue. Distant metastasisoccurs late during the genetic evolution of pancreatic cancer. Nature, 467(7319):1114–7, 2010. → pages 3, 19[137] L. R. Yates, M. Gerstung, S. Knappskog, C. Desmedt, G. Gundem, P. VanLoo, T. Aas, L. B. Alexandrov, D. Larsimont, H. Davies, Y. Li, Y. S. Ju,M. Ramakrishna, H. K. Haugland, P. K. Lilleng, S. Nik-Zainal,S. McLaren, A. Butler, S. Martin, D. Glodzik, A. Menzies, K. Raine,J. Hinton, D. Jones, L. J. Mudie, B. Jiang, D. Vincent, A. Greene-Colozzi,P.-Y. Adnet, A. Fatima, M. Maetens, M. Ignatiadis, M. R. Stratton,C. Sotiriou, A. L. Richardson, P. E. Lønning, D. C. Wedge, and P. J.Campbell. Subclonal diversification of primary breast cancer revealed bymultiregion sequencing. Nature medicine, 21(7):751–759, 2015. → page 3[138] C. Yu, J. Yu, X. Yao, W. K. K. Wu, Y. Lu, S. Tang, X. Li, L. Bao, X. Li,Y. Hou, R. Wu, M. Jian, R. Chen, F. Zhang, L. Xu, F. Fan, J. He, Q. Liang,H. Wang, X. Hu, M. He, X. Zhang, H. Zheng, Q. Li, H. Wu, Y. Chen,147X. Yang, S. Zhu, X. Xu, H. Yang, J. Wang, X. Zhang, J. J. Y. Sung, Y. Li,and J. Wang. Discovery of biclonal origin and a novel oncogene SLC12A5in colon cancer by single-cell sequencing. Cell Research, 24(6):701–712,2014. → page 122[139] K. Yuan, T. Sakoparnig, F. Markowetz, and N. Beerenwinkel.Bitphylogeny: a probabilistic framework for reconstructing intra-tumorphylogenies. Genome biology, 16(1):36, 2015. → page 19[140] H. Zafar, Y. Wang, L. Nakhleh, N. Navin, and K. Chen. Monovar:single-nucleotide variant detection in single cells. Nature Methods, 13:505,apr 2016. → pages 15, 58[141] H. Zahn, A. Steif, E. Laks, P. Eirew, M. VanInsberghe, S. P. Shah,S. Aparicio, and C. L. Hansen. Scalable whole-genome single-cell librarypreparation without preamplification. Nature Methods, 14(2):167–173,Feb. 2017. → pages 10, 11, 56, 57, 60, 84[142] A. W. Zhang, A. McPherson, K. Milne, D. R. Kroeger, P. T. Hamilton,A. Miranda, T. Funnell, N. Little, C. P. de Souza, S. Laan, et al. Interfacesof malignant and immunologic clonal dynamics in ovarian cancer. Cell,173(7):1755–1769, 2018. → page 2[143] J. Zhang, J. Fujimoto, J. Zhang, D. C. Wedge, X. Song, J. Zhang, S. Seth,C.-W. Chow, Y. Cao, C. Gumbs, K. A. Gold, N. Kalhor, L. Little,H. Mahadeshwar, C. Moran, A. Protopopov, H. Sun, J. Tang, X. Wu, Y. Ye,W. N. William, J. J. Lee, J. V. Heymach, W. K. Hong, S. Swisher, I. I.Wistuba, and P. A. Futreal. Intratumor heterogeneity in localized lungadenocarcinomas delineated by multiregion sequencing. Science (NewYork, N.Y.), 346(6206):256–9, 2014. → page 3[144] Y. Zhou, A. M. Johansen, and J. A. Aston. Toward automatic modelcomparison: An adaptive sequential monte carlo approach. Journal ofComputational and Graphical Statistics, 25(3):701–726, 2016. → page 64[145] C. Zong, S. Lu, A. R. Chapman, and X. S. Xie. Genome-wide detection ofsingle-nucleotide and copy-number variations of a single human cell.Science, 338(6114):1622–1626, 2012. → page 56148
Thesis/Dissertation
eng
Computer Science
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Somatic mutation analysis for the study of clonal evolution in cancer
Text