ipcoal

What it simulates

ipcoal simulates gene genealogies and DNA variants along a species tree under the multispecies coalescent. It builds on msprime for the simulation and adds a species-tree interface, per-locus genealogy tables and tools for phylogenetic analysis. You give it a species tree with branch lengths in generations, an effective population size, the number of samples per tip, the number of loci and sites, and mutation and recombination rates.

It returns the local gene tree for every stretch of sequence between recombination breakpoints, plus the SNPs that fell on them. The key assumption is that gene trees can differ from the species tree through incomplete lineage sorting: two lineages that fail to coalesce inside a short branch may coalesce deeper, in a different order. The recipe's helper script uses one Ne for every branch. Loci are unlinked from each other and recombine only inside a locus.

When to use it

Use ipcoal to ask how often gene trees disagree with the species tree for a given Ne and branch length. It also generates test data for species-tree and gene-tree inference methods, and for checking a phylogenetic pipeline against a known truth.

• Use msprime when you have one population or a custom demography and want a VCF with full control over the model.
• Use stdpopsim for published models of a real species.

Run it in your browser

ipcoal has no browser version on biosimhub. Use the install route below.

Install it yourself

The conda command resolves for ipcoal 0.5.0:

micromamba create -n ipcoal -c conda-forge ipcoal=0.5.0

The ipcoal docs also list pip install ipcoal. A minimal run in Python, without the recipe's helper script:

import ipcoal, toytree

tree = toytree.tree("((A:200000,B:200000):200000,(C:200000,D:200000):200000);")
model = ipcoal.Model(tree, Ne=100_000, nsamples=2, mut=1e-8, recomb=1e-9,
                     seed_trees=1, seed_mutations=1)
model.sim_loci(nloci=3, nsites=1000)
print(model.df[["locus", "start", "end", "nsnps"]])

In the recipe's 0.5.0 image, locus 0 split into 5 genealogy segments with 3, 11, 9, 3 and 0 SNPs. The table model.df has one row per segment.

Reading the output

The recipe's helper script run_ipcoal.py reads your Newick tree with toytree, builds ipcoal.Model with the same seed for trees and mutations, calls sim_loci, writes the VCF with model.write_vcf, and saves model.df as CSV. The example tree is ((A:200000,B:200000):200000,(C:200000,D:200000):200000);. It has four tips, and each tip is sampled twice, so the data have 8 sequences named A_0, A_1, B_0 and so on.

• ipcoal_log.txt: "done 19 genealogy segments, 274 SNPs". Five loci of 2000 sites each, cut by recombination into 19 segments, with 274 SNPs.
• genealogies.csv: columns locus,start,end,nbps,nsnps,tidx,genealogy. Locus 0 has two rows. Positions 0 to 115 (115 bp) carry 1 SNP, and positions 115 to 2000 (1885 bp) carry 40 SNPs, each with its own Newick tree. Locus 1 splits at 521 and 844.
• gene_trees.nwk: the same 19 trees, one per line. Branch lengths are in generations. The first tree joins B_0 and B_1 at 202,832 generations and C_0 and C_1 at 29,804.
• ipcoal.vcf: contigs 0 to 4 of length 2000, with 8 sample columns (A_0, A_1, B_0, B_1, C_0, C_1, D_0, D_1).
• summary.txt: seed, the tree, 5 loci, 2000 sites per locus, 19 segments and 274 SNPs.

The species tree groups A with B and C with D. The first gene tree has the B pair outermost and the second has the C pair outermost, so both disagree with it. Gene trees that differ from the species tree are incomplete lineage sorting, the effect ipcoal is built to study.

Your own runs

Each run changes one value from the example (Ne 100,000, 2 samples per tip, 5 loci, seed 42).

• Smaller Ne. --ne 10000 gave 5 genealogy segments and 141 SNPs, down from 19 and 274. With a tenth of the Ne there are fewer recombination breakpoints and fewer SNPs, because both depend on coalescent time, which scales with Ne.
• No recombination. --recomb 0 gave 5 segments (one per locus) and 309 SNPs.
• More loci. --nloci 20 gave 100 segments and 1010 SNPs. Four times the loci gave about four times the SNPs.
• More samples per tip. --nsamples 4 gave 24 segments and 341 SNPs, with 16 sequences in the VCF.

Pitfalls

• Branch lengths are generations. The tree you give is in generations. Branch lengths in substitutions per site would make every split a fraction of a generation old.
• Mutation rate is per site per generation. The default 1e-8 gives about 2.7 SNPs per 100 sites for the example (274 over 10,000 sites). Short loci need a larger rate or more sites to have any SNPs.
• Loci are independent. Recombination acts within a locus. For linked sequence, use one locus with more sites.
• Same seed for trees and mutations. The recipe's helper script uses one seed for both. Change --seed to get a new replicate.
• Tip names set the sample names. Sequences are named <tip>_<n>, so a tip called "A" is sampled as A_0, A_1.

Further reading

• Documentation: https://eaton-lab.org/ipcoal/
• Source: https://github.com/eaton-lab/ipcoal
• Paper: McKenzie PF, Eaton DAR (2020) ipcoal: an interactive Python package for simulating and analyzing genealogies and sequences on a species tree or network. Bioinformatics 36:4193-4196. https://doi.org/10.1093/bioinformatics/btaa486