search

daijiang / phyr

Functions for phylogenetic analyses
https://daijiang.github.io/phyr/
GNU General Public License v3.0
30 stars 10 forks source link

pcd function not working with ultrametric tree #49

Closed graciellehigino closed 4 years ago

graciellehigino commented 4 years ago

Hello!

Thank you for maintaining this package, it's really useful. I am having a problem when I try to run the function pcd with an ultrametric tree. I need to calculate PCD for two trees using the same community matrix, and it works for one of them, but not for the other. I already checked if all species of the tree tips are in the community matrix and vice-versa. When I tried to run line by line and setting the cpp argument as FALSE, it runs, but the result is full of NA's (and it was not supposed to be like that). Am I missing something?

Thank you again!

library(ape)
library(phyr)

hp_mat <- read.csv("https://gist.githubusercontent.com/kguidonimartins/9955c6d30733a9b4a196c43289f467da/raw/3a11137277a2fdfb3b414176409067dbaab30d6f/hp_mat.csv", sep=";")

tree_fleas <- read.tree("https://gist.githubusercontent.com/kguidonimartins/9955c6d30733a9b4a196c43289f467da/raw/3a11137277a2fdfb3b414176409067dbaab30d6f/FleaTree.tre")

tree_hosts <- read.tree("https://gist.githubusercontent.com/kguidonimartins/9955c6d30733a9b4a196c43289f467da/raw/3a11137277a2fdfb3b414176409067dbaab30d6f/HostTree.tre")

pcd_fleas <- pcd(hp_mat, tree_fleas)
#> Dropping tips from the phylogeny that are not in the community data
#> Dropping species from the community data that are not in the phylogeny
#> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

pcd_hosts <- pcd(hp_mat, tree_hosts)
#> Dropping tips from the phylogeny that are not in the community data
#> Dropping species from the community data that are not in the phylogeny
#> 1
#> Error in pcd2_loop(SSii, nsr, SCii, as.matrix(comm), V, nsp_pool, verbose): Mat::submat(): indices out of bounds or incorrectly used
daijiang commented 4 years ago

Thanks! Github changed its notification center, which I find is more annoying. Thus I did not see your issue till now.

The problem here is that one of your site, after removing species that are not in the phylogeny, has no species anymore.

hp_mat2 = hp_mat[, colnames(hp_mat) %in% tree_hosts$tip.label]
which(rowSums(hp_mat2) == 0) # 22

If you remove this site, it should work. I have updated the function in case this happens later. Thanks!

graciellehigino commented 4 years ago

Yaass thank you! It worked perfectly! [=