ropensci / phylocomr

Phylocom R interface
https://docs.ropensci.org/phylocomr
Other
15 stars 7 forks source link

fixed node ages in ph_bladj are apparently not taken into account #19

Closed LunaSare closed 7 years ago

LunaSare commented 7 years ago
Session Info ``` Session info ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- setting value version R version 3.4.1 (2017-06-30) system x86_64, darwin15.6.0 ui AQUA language (EN) collate en_US.UTF-8 tz America/New_York date 2017-10-04 Packages ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- package * version date source ape * 4.1 2017-02-14 CRAN (R 3.4.0) base * 3.4.1 2017-07-07 local coda 0.19-1 2016-12-08 CRAN (R 3.4.0) compiler 3.4.1 2017-07-07 local curl 2.8.1 2017-07-21 CRAN (R 3.4.1) datasets * 3.4.1 2017-07-07 local deSolve 1.20 2017-07-14 CRAN (R 3.4.1) devtools 1.13.3 2017-08-02 CRAN (R 3.4.1) digest 0.6.12 2017-01-27 CRAN (R 3.4.0) geiger * 2.0.6 2015-09-07 CRAN (R 3.4.0) git2r 0.19.0 2017-07-19 CRAN (R 3.4.1) graphics * 3.4.1 2017-07-07 local grDevices * 3.4.1 2017-07-07 local grid 3.4.1 2017-07-07 local httr 1.3.1 2017-08-20 cran (@1.3.1) lattice 0.20-35 2017-03-25 CRAN (R 3.4.1) MASS 7.3-47 2017-02-26 CRAN (R 3.4.1) memoise 1.1.0 2017-04-21 CRAN (R 3.4.0) methods * 3.4.1 2017-07-07 local mvtnorm 1.0-6 2017-03-02 CRAN (R 3.4.0) nlme 3.1-131 2017-02-06 CRAN (R 3.4.1) parallel 3.4.1 2017-07-07 local phylocomr * 0.1.0 2017-09-26 Github (ropensci/phylocomr@3d2eafa) R6 2.2.2 2017-06-17 CRAN (R 3.4.0) Rcpp 0.12.12 2017-07-15 CRAN (R 3.4.1) stats * 3.4.1 2017-07-07 local subplex 1.4-1 2017-07-18 CRAN (R 3.4.1) sys 1.4 2017-06-24 CRAN (R 3.4.1) utils * 3.4.1 2017-07-07 local withr 2.0.0 2017-07-28 CRAN (R 3.4.1) ```

Using the following commands, with three diferent values of fixed internal node ages, we get the same output trees:

devtools::install_github("ropensci/phylocomr")
library(phylocomr)
library(geiger)
Loading required package: ape
ages_file <- system.file("examples/ages", package = "phylocomr")
phylo_file <- system.file("examples/phylo_bladj", package = "phylocomr")
# from data.frame
ages_df <- data.frame(
   a = c('malpighiales','salicaceae','fabaceae','rosales','oleaceae',
         'gentianales','apocynaceae','rubiaceae'),
   b = c(81,20,56,76,47,71,18,56)
 )
phylo_str <- readLines(phylo_file)
res <- phylocomr::ph_bladj(ages = ages_df, phylo = phylo_str)
plot(read.tree(text = res))
nodelabels()
ape::axisPhylo()
ages_df_me <- data.frame(
   a = c('malpighiales','salicaceae','fabaceae','rosales','oleaceae',
         'gentianales','apocynaceae','rubiaceae'),
   b = c(81,120,156,176,147,171,118,156)
 )
res2 <- phylocomr::ph_bladj(ages = ages_df_me, phylo = phylo_str) # 
plot(read.tree(text = res2))
nodelabels()
ape::axisPhylo()
res
# [1] "((((((lomatium_concinnum:20.250000,campanula_vandesii:20.250000):20.250000,(((veronica_candidissima:10.125000,penstemon_paniculatus:10.125000)plantaginaceae:10.125000,justicia_oblonga:20.250000):10.125000,marsdenia_gilgiana:30.375000):10.125000):10.125000,epacris_alba-compacta:50.625000)ericales_to_asterales:10.125000,((daphne_anhuiensis:20.250000,syzygium_cumini:20.250000)malvids:20.250000,ditaxis_clariana:40.500000):20.250000):10.125000,thalictrum_setulosum:70.875000)eudicots:10.125000,((dendrocalamus_giganteus:27.000000,guzmania_densiflora:27.000000)poales:27.000000,warczewiczella_digitata:54.000000):27.000000)malpighiales:1.000000;\n"
# attr(,"ages_file")
# [1] "/var/folders/ss/2tpkp325521_kfgn59g44vd80000gn/T//RtmpaZxQEc/ages"
# attr(,"phylo_file")
# [1] "/var/folders/ss/2tpkp325521_kfgn59g44vd80000gn/T//RtmpaZxQEc/phylo_421b1322db5"
res2
# [1] "((((((lomatium_concinnum:20.250000,campanula_vandesii:20.250000):20.250000,(((veronica_candidissima:10.125000,penstemon_paniculatus:10.125000)plantaginaceae:10.125000,justicia_oblonga:20.250000):10.125000,marsdenia_gilgiana:30.375000):10.125000):10.125000,epacris_alba-compacta:50.625000)ericales_to_asterales:10.125000,((daphne_anhuiensis:20.250000,syzygium_cumini:20.250000)malvids:20.250000,ditaxis_clariana:40.500000):20.250000):10.125000,thalictrum_setulosum:70.875000)eudicots:10.125000,((dendrocalamus_giganteus:27.000000,guzmania_densiflora:27.000000)poales:27.000000,warczewiczella_digitata:54.000000):27.000000)malpighiales:1.000000;\n"
# attr(,"ages_file")
# [1] "/var/folders/ss/2tpkp325521_kfgn59g44vd80000gn/T//RtmpaZxQEc/ages"
# attr(,"phylo_file")
# [1] "/var/folders/ss/2tpkp325521_kfgn59g44vd80000gn/T//RtmpaZxQEc/phylo_421b3ff2e868"
ages_df_me <- data.frame(
   a = c('malpighiales','salicaceae','fabaceae','rosales','oleaceae',
         'gentianales','apocynaceae','rubiaceae'),
   b = c(81,12.0,15.6,17.6,14.7,17.1,11.8,15.6)
 )
res3 <- phylocomr::ph_bladj(ages = ages_df_me, phylo = phylo_str) # 
plot(read.tree(text = res3))
nodelabels()
ape::axisPhylo()
res3
# [1] "((((((lomatium_concinnum:20.250000,campanula_vandesii:20.250000):20.250000,(((veronica_candidissima:10.125000,penstemon_paniculatus:10.125000)plantaginaceae:10.125000,justicia_oblonga:20.250000):10.125000,marsdenia_gilgiana:30.375000):10.125000):10.125000,epacris_alba-compacta:50.625000)ericales_to_asterales:10.125000,((daphne_anhuiensis:20.250000,syzygium_cumini:20.250000)malvids:20.250000,ditaxis_clariana:40.500000):20.250000):10.125000,thalictrum_setulosum:70.875000)eudicots:10.125000,((dendrocalamus_giganteus:27.000000,guzmania_densiflora:27.000000)poales:27.000000,warczewiczella_digitata:54.000000):27.000000)malpighiales:1.000000;\n"
# attr(,"ages_file")
# [1] "/var/folders/ss/2tpkp325521_kfgn59g44vd80000gn/T//RtmpaZxQEc/ages"
# attr(,"phylo_file")
# [1] "/var/folders/ss/2tpkp325521_kfgn59g44vd80000gn/T//RtmpaZxQEc/phylo_421b93ea06d"
sckott commented 7 years ago

thanks @LunaSare for the issue!

will have a look!

sckott commented 7 years ago

@LunaSare

So I think the example I had for bladj was just not a good example as some of the node names in the data.frame in the example were not even in the tree in the example. I'm redoing the documentation for ph_bladj now, but here's example that should work:

library(phylocomr)
phylo_file <- system.file("examples/phylo_bladj", package = "phylocomr")
phylo_str <- readLines(phylo_file)
df <- data.frame(
    a = c('malpighiales','eudicots','ericales_to_asterales','plantaginaceae', 'malvids', 'poales'),
    b = c(81,20,56,76,47,71)
)
res <- ph_bladj(ages = df, phylo = phylo_str)
attributes(res) <- NULL

df2 <- data.frame(
    a = c('malpighiales','eudicots','ericales_to_asterales','plantaginaceae', 'malvids', 'poales'),
    b = c(91,21,54,75,49,73)
)
res1 <- ph_bladj(ages = df2, phylo = phylo_str)
attributes(res1) <- NULL
identical(res, res1)
#> FALSE
res
#> [1] "((((((lomatium_concinnum:5.714286,campanula_vandesii:5.714286):5.714286,(((veronica_candidissima:2.857143,penstemon_paniculatus:2.857143)plantaginaceae:2.857143,justicia_oblonga:5.714286):2.857142,marsdenia_gilgiana:8.571428):2.857143):2.857143,epacris_alba-compacta:14.285715)ericales_to_asterales:2.857141,((daphne_anhuiensis:5.714285,syzygium_cumini:5.714285)malvids:5.714285,ditaxis_clariana:11.428571):5.714286):2.857143,thalictrum_setulosum:20.000000)eudicots:61.000000,((dendrocalamus_giganteus:71.000000,guzmania_densiflora:71.000000)poales:5.000000,warczewiczella_digitata:76.000000):5.000000)malpighiales:1.000000;\n"

res1
#> [1] "((((((lomatium_concinnum:6.000000,campanula_vandesii:6.000000):6.000000,(((veronica_candidissima:3.000000,penstemon_paniculatus:3.000000)plantaginaceae:3.000000,justicia_oblonga:6.000000):3.000000,marsdenia_gilgiana:9.000000):3.000000):3.000000,epacris_alba-compacta:15.000000)ericales_to_asterales:3.000000,((daphne_anhuiensis:6.000000,syzygium_cumini:6.000000)malvids:6.000000,ditaxis_clariana:12.000000):6.000000):3.000000,thalictrum_setulosum:21.000000)eudicots:70.000000,((dendrocalamus_giganteus:73.000000,guzmania_densiflora:73.000000)poales:9.000000,warczewiczella_digitata:82.000000):9.000000)malpighiales:1.000000;\n"
sckott commented 7 years ago

@LunaSare i updated the documentation, can reinstall to get it, like devtools::install_github("ropensci/phylocomr")

LunaSare commented 7 years ago

thanks @sckott! the examples seem to be working fine now. I also tried just changing internal nodes and it gives different outputs now:

devtools::install_github("ropensci/phylocomr") library(phylocomr) phylo_file <- system.file("examples/phylo_bladj", package = "phylocomr") phylo_str <- readLines(phylo_file) df <- data.frame( a = c('malpighiales','eudicots','ericales_to_asterales','plantaginaceae', 'malvids', 'poales'), b = c(81,20,56,76,47,71) ) res <- ph_bladj(ages = df, phylo = phylo_str) attributes(res) <- NULL df3 <- data.frame( a = c('malpighiales','eudicots','ericales_to_asterales','plantaginaceae', 'malvids', 'poales'), b = c(81,21,54,75,49,73) ) res3 <- ph_bladj(ages = df3, phylo = phylo_str) attributes(res3) <- NULL identical(res,res3) [1] FALSE

sckott commented 7 years ago

great! glad it works.