Estimating Mutation Parameters and Population History Simultaneously from Temporally-Spaced Genome Data Arman Bilge, Tanja Stadler, Matthew Kearse, and Alexei J. Drummond email: abil933@aucklanduni.ac.nz Estimating Mutation Parameters and Population History Simultaneously from Temporally-Spaced Genome Data Arman Bilge, Tanja Stadler, Matthew Kearse, and Alexei J. Drummond email: abil933@aucklanduni.ac.nz Motivation and Primary Challenges ▷Very feasible to sequence entire genomes ▷More recently, even possible to recover ancient genomes ▷Temporally-spaced genome data ▷Opportunity to do inference previously only possible for fast-evolving organisms (e.g., viruses) ▷mutation rate ▷population size through time But. . . ▷Diﬃcult to phase diploid genomes ▷Recombination ▷Low coverage and sequencing error, especially for ancient genomes ▷Cannot use existing Bayesian phylogenetic methods Sequencing Depth and Error ▷Assume that all individuals have about same number of SNPs ▷Average sequencing depth of sample is correlated with observed SNPs ▷Our dataset approaches complete SNP recovery at 22x coverage ▷Ancient genomes are sequenced at lower depth and thus missing many SNPs ▷Leads to systematic bias in estimates 4 8 12 16 20 24 28 4 4.5 5 5.5 6 6.5 7 ·106 Average Depth Recovered SNPs Overview of Methodology ▷Want to estimate mutation and population parameters θ from pileup data ▷pileup data is unsummarized, aligned reads for each sampled individual ▷To compute posterior need to marginalize over individual’s genotypes ▷Computationally intractable so use importance sampling ▷The importance distribution assumes independence of individual’s genotypes P (θ | D) = X G P (θ | G) P (G | D) = lim n→∞ n X i=1 P  θ | G (i) P G (i) | D  ˆP G (i) | D , G (i) ∼ˆP (· | D) ▷Finally, use standard MCMC to sample parameters for a given genotype P (θ | G) ∝P (G | θ) P (θ) Sampling an Individual’s Genotype ▷Want the genotype g1, g2 ∈{A, C, G, T} of individual at a position in its genome ▷Data is the observed base calls at this position with their Phred quality scores D = (bq : b ∈{A, C, G, T} , q ∈N) ▷Number of base calls |D| is the sequencing depth at this position ▷Sample genotype from the posterior distribution P (g1, g2 | D) = P (D | g1, g2) P (g1) P (g2) P (D) ▷To compute P (D | g1, g2) assume base calls are multinomially distributed with probabilities P (bq | g1, g2) = 1 2P (bq | g1) + 1 2P (bq | g2) ▷Using the deﬁnition of a Phred quality score (assuming equal error rates for all bases) P (gi | bq) = ( 1 −10 −q 10 ifb = gi 1 310 −q 10 ifb ̸= gi we have P (bq | gi) = P (gi | bq) P (bq) P (gi) ▷Use empirical estimates for P (gi) and P (bq) Genotype Probability ▷Assumes that sites are unlinked; i.e., phylogenetically independent ▷Models recombination with no dependence on correct phasing ▷Assumes a site is biallelic; i.e., has only two possible nucleotide states ▷Often true in practice ▷Can be handled rigorously with ascertainment bias correction ▷Want the probability of all the individuals’ genotypes by marginalizing over all phylogenies P (G | θ) = Z T P (G | T, θ) P (T | θ) dT ▷Integral is computed numerically using similar technique to SNAPP [2] ▷Divide time into intervals using sampling times and population change times ▷Each interval i can be described by a linear system of diﬀerential equations Qi ▷Solve each system by taking the matrix exponential exp Qi [1] ▷Can be done eﬃciently by caching and reusing matrix exponentials Simulation Study ▷8 diploid taxa, including 4 ancient individuals up to 50k years old ▷Mutation rate µ = 10−9 s/s/yr ▷HKY model with κ = 5 ▷Constant size population with Ne = 3 × 106 ▷Simulated 50 datasets of 104 total sites ▷Attempted to infer parameters ▷True values always within 95% HPD ▷Mean ˆµ within ±1.8 × 10−10 of true µ G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G 0e+00 1e−09 2e−09 3e−09 0 10 20 30 40 50 replicates Mean Summary ▷Fully Bayesian inference of mut