Maya A. Lewinsohn1,2, Trevor Bedford1,2,3, Nicola F. Müller2, Alison F. Feder1
1Department of Genome Sciences, University of Washington, Seattle, WA, USA; 2Vaccine and Infectious Disease Division, Fred Hutchinson Cancer Research Center, Seattle, WA, USA; 3Howard Hughes Medical Institute, Seattle, WA, USA
Spatial properties of tumor growth have profound implications for cancer progression, therapeutic resistance and metastasis. Yet, how spatial position governs tumor cell division remains difficult to evaluate in clinical tumors. Here, we demonstrate that elevated cellular growth rates on the tumor periphery leave characteristic patterns in the genomes of cells sampled from different parts of a tumor, which become evident when they are used to construct a tumor phylogenetic tree. Namely, rapidly-dividing peripheral lineages branch more extensively and acquire more mutations than slower-dividing lineages in the tumor center. We develop a Bayesian state-dependent evolutionary phylodynamic model (SDevo) that quantifies these patterns to infer the differential cell division rates between peripheral and central cells jointly from the branching and mutational patterns of single-time point, multi-region sequencing data. We validate this approach on simulated tumors by demonstrating its ability to accurately infer spatially-varying birth rates under a range of growth conditions and sampling strategies. We then show that SDevo outperforms state-of-the-art, non-cancer multi-state phylodynamic methods which ignore differential mutational acquisition. Finally, we apply SDevo to multi-region sequencing data from clinical hepatocellular carcinomas and find evidence that cells on the tumor edge divide 3-6x faster than those in the center. As multi-region and single-cell sequencing increase in resolution and availability, we anticipate that SDevo will be useful in interrogating spatial restrictions on tumor growth and could be extended to model non-spatial factors that influence tumor progression, including hypoxia and immune infiltration.
Scripts to generate simulated tumors under boundary-driven and unrestricted growth in a 2D lattice can be found in simulated_data/spatial_tumor_simulation.ipynb.
Simulated tumors are recorded .csv files recording all cells in tumor simulation. For pushing simulation, cell positions in lattice are
recorded in locs.csv at each time slice.
scripts/reconstruct_simulated_trees.R contains code to convert simulated cells records into S4 tree objects containing edge/center state information.
Install local R package tumortree for necessary functions.
Simulation data is generated from https://github.com/federlab/PhysiCellTrees and should be put phyiscell/simulation_data for downstream analyses.
Script to generate input XML files from simulated Eden simulation can be found here: scripts/set_up_simulated_tumors_xmls.R
Script to generate input XML files from PhysiCell outputs can be found here: scripts/write_state_clocks_xml_from_physicell.R.
run_sdevo_simulation_study.sh provides commands for running SDevo on all XML files using a Slurm workload manager. Essentially, run_beast_to_ess.sh is called until each MCMC chain reaches the desired ESS. check_ess_mcmc.Rscript is used here to check the current ESS of a given output log and write summary statistics of the posterior estimates for each run. compile_posterior_summaries.sh can be run at the end to generate TSV combining all summary statistics.
Make XML files with strict clock by running create_strict_clock.R, where input is state clocks XML file generated above.
Scripts to generate input DNA sequences HCC sequence data can be found in scripts/process_li_data.R and li-application/. Maximum likelihood trees are solved by Fasttree and Augur, see run_nextstrain_divergence_trees.sh.
Local R package tumortree is needed for most figures.