Find out unit of phangorn::simSeq

[richelbilderbeek]

This is what the doc says:

rate mutation rate or scaler for the edge length, a numerical value greater than zero.

I would have liked to read:

rate mutation rate, nucleotide substitution rate per unit of time

This code generates an alignment with maximum information (mutation rate = 1K, time = 15, DNA alignment length = 15K), but it appears all garbled:

out <- flechette::run(
  flechette::create_params(
    speciation_initiation_rate = 0.1,
    speciation_completion_rate = 1.0,
    extinction_rate = 0.01,
    crown_age = 15,
    crown_age_sigma = 0.001,
    sampling_method = "youngest",
    mutation_rate = 1000,
    sequence_length = 15000,
    mcmc_length = 10000,
    tree_sim_rng_seed = 42,
    alignment_rng_seed = 42,
    beast2_rnd_seed = 42
  ),
  verbose = TRUE
)
names(out)
names(out$estimates)
image(out$alignment)

filename <- tempfile()
phangorn::write.phyDat(out$alignment, file = filename, format = "fasta")

tree <- phangorn::optim.parsimony(
  tree = phangorn::random.addition(phangorn::as.phyDat(x = out$alignment)),
  data = phangorn::as.phyDat(x = out$alignment)
)
ape::plot.phylo(tree)
ape::plot.phylo(out$species_tree)

I don't mind considering being 'just human' and unable to see the pattern that is there, but I want to be sure.

[richelbilderbeek]

Bug was at my side: I uses a mutation rate of 1000 bp per million years per lineage. As I use 15 million years, one expects 15 (some silent) mutations per base pair; also called 'a meaningless jumble of nucleotides'.