thijsjanzen / treestats

https://thijsjanzen.github.io/treestats/
14 stars 2 forks source link

Apply treestats::crown_age() to non-ultrametric tree gives inconsistent results, and more #20

Open EvoLandEco opened 1 year ago

EvoLandEco commented 1 year ago

DATA Original data of the tree brv_178_sm_appendixs5.xls, from A phylogeny of Cenozoic macroperforate planktonic foraminifera from fossil data

A tree export to txt (github does not support uploading a .phy) I converted from sheet aMb (fully bifurcating morphospecies phylogeny) using their paleoPhylo code tree_morpho_aMb.txt

This is the code to convert the sheet data to phylo tree, the functions required are in paleoPhylo.txt

tree_morpho <-
  with(
    tree_morpho_table,
    as.paleoPhylo(`Species code`, `Ancestor code`, `Start date`, `End date`, label =
                    `Species name`)
  )
tree_morpho <- buildApe(createBifurcate(tree_morpho))

This is the complete tree: image

EXAMPLE

tree_morpho <- ape::read.tree("c:/tree_morpho_aMb.txt")

# beautier does not support non-ultrametric tree
beautier::get_crown_age(tree_morpho)
#> Error in beautier::get_crown_age(tree_morpho): 'phylogeny' must be ultrametric

# if we drop extinct species, it works
beautier::get_crown_age(geiger::drop.extinct(tree_morpho))
#> [1] 65

# apply treestats::crown_age() to tree without extinct lineages gives the same result
treestats::crown_age(geiger::drop.extinct(tree_morpho))
#> [1] 65

# apply treestats::crown_age() to the complete tree gives inconsistent result
treestats::crown_age(tree_morpho)
#> [1] 69.2

# test again by a tree generated by DDD
tree_ddd <- DDD::dd_sim(pars = c(0.6,0.2, 100), age = 50)

# both beautier and treestats work as expected
beautier::get_crown_age(tree_ddd$tas)
#> Error in beautier::get_crown_age(tree_ddd$tas): 'phylogeny' must be ultrametric
beautier::get_crown_age(tree_ddd$tes)
#> [1] 50

treestats::crown_age(tree_ddd$tes)
#> [1] 50
treestats::crown_age(tree_ddd$tas)
#> [1] 50

Created on 2023-08-15 with reprex v2.0.2

MORE Also interesting that before I exported the tree to newick format and read it back again, applying crown_age() to the just-converted tree gave a different result

> treestats::crown_age(tree_morpho_bak)
[1] 9.9

By furthur inspecting the results, I found that the two trees tree_morpho (the one exported to newick and read back) and tree_morpho_bak (originally converted by paleoPhylo) also differed in the following stats:

all_stats_a <- treestats::calc_all_stats(tree_morpho)
all_stats_b <- treestats::calc_all_stats(tree_morpho_bak)
setdiff(all_stats_a, all_stats_b)

$gamma
[1] 414.1048

$beta
[1] -1.390032

$crown_age
[1] 69.2

$pigot_rho
[1] -0.2383187

$nltt_base
[1] 0.9911671

$phylogenetic_div
[1] 684.47

$laplac_spectrum_a
[1] 0.1448609

$var_depth
[1] 46.43248

$mpd
[1] 68.40613

$psv
[1] 34.20306

$vpd
[1] 891.51

$eigenvector
[1] 0.7071068

But if extinct linneages were dropped, the two trees only differed in four stats:

all_stats_a2 <- treestats::calc_all_stats(geiger::drop.extinct(tree_morpho))
all_stats_b2 <- treestats::calc_all_stats(geiger::drop.extinct(tree_morpho_bak))
setdiff(all_stats_a2, all_stats_b2)

$laplac_spectrum_a
[1] 0.3689392

$laplac_spectrum_p
[1] 2.071078

$psv
[1] 48.8603

$vpd
[1] 1974.94
EvoLandEco commented 1 year ago

It seems that paleoPhylo generated a mal-formatted tree, but I think it is worthy to check if it could be a rare case that causes error

thijsjanzen commented 1 year ago

It seems that paleoPhylo generated a mal-formatted tree, but I think it is worthy to check if it could be a rare case that causes error

I'm not sure how to replicate this error, or what now caused the documented issue - I also can't really replicate the tree plotted here, what package did you use to read the weird newick string?

EvoLandEco commented 1 year ago

Hi @thijsjanzen, did you try the code in paleoPhylo.txt? This is the code the authors provided to convert their data (table stored in xlsx) into a phylo object.

thijsjanzen commented 1 year ago

Right, I'm not sure what happens during your newick conversion-reread, so I can't check that.

However, I do notice that treestats::crown_age provides a weird result when applied to the full tree (9.9, instead of 65). Dunno yet was is going on, interesting!