ContMap throws error for some trees

[voidikova]

Hello, I am using the contMap function but it reports this error for some trees:

Error in if (x <= trans[1]) state <- names(trans)[1] else if (x >= trans[length(trans)]) state <- names(trans)[length(trans)] else { : 
  missing value where TRUE/FALSE needed

Below is the example for two trees. They differ only in one species (Lunaria rediviva is included in rawBadTree and missing in the rawGoodTree).

library(phytools)

testContMap <- function(data) {
    tree <- read.tree(text=data)
    traits <- rnorm(length(tree$tip.label))
    names(traits) <- tree$tip.label
    contMap(tree, traits)
}

testNodes <- function(data) {
    tree <- read.tree(text=data)
    btree <- multi2di(tree)
    head(matchNodes(tree, btree, method = "distances", quiet = TRUE), n=2)
}

rawBadTree = "((Juncus_conglomeratus:73.524994,((Carex_ovalis:26.6875,((Carex_vulpina:5.71875,Carex_otrubae:5.71875)Carex_vulpina_agg.:17.15625,Carex_diandra:22.875)N1205:3.8125)N1190:6.8125,(Carex_tomentosa:12.300001,(Carex_melanostachya:10.2,(Carex_pseudocyperus:9.439999,(Carex_distans:8.1,(Carex_pendula:5.3,Carex_demissa:5.3)N1379:2.8)N1370:1.339999)N1326:0.76)N1313:2.100001)N1256:21.199999)N1185:40.024992)N1063:74.27501,(Chelidonium_majus:135,((((Hypericum_montanum:5,Hypericum_hirsutum:5)N2928:98.099997,Astragalus_glycyphyllos:103.099998)N2698:5.900001,((Descurainia_sophia:27.299999,((Lepidium_campestre:13.2,Coronopus_squamatus:13.2)Lepidieae:12.373334,(Camelina_sativa:6.6,Camelina_microcarpa:6.6,Camelina_alyssum:6.6)Camelina:18.973334)N4741:1.726666)Lineage_I_Franzke_et_al._2011:5,Bunias_orientalis:32.300001,Lunaria_rediviva:32.299998,((Berteroa_incana:19.3,Alyssum_alyssoides:19.300001)Alysseae:12.325,Arabis_hirsuta:31.625,(Isatis_tinctoria:29.600002,Diplotaxis_tenuifolia:29.6)Lineage_II_Franzke_et_al._2011:2.024999)Expanded_Lineage_II_Franzke_et_al._2011:0.675)Core_Brassicaceae:76.700001)Rosidae:9.199999,((((Dianthus_deltoides:4.53,Dianthus_gratianopolitanus:4.53)Dianthus:27.04,Holosteum_umbellatum:31.569999)N5370:27.755001,(Bassia_scoparia:51.730001,Atriplex_rosea:51.73)N5582:7.594999)N5289:52.175002,Adenostyles_alliariae:111.500003)N5136:6.7)N2487:16.800009)Eudicots:12.800003)N398;"

rawGoodTree = "((Juncus_conglomeratus:73.524994,((Carex_ovalis:26.6875,((Carex_vulpina:5.71875,Carex_otrubae:5.71875)Carex_vulpina_agg.:17.15625,Carex_diandra:22.875)N1205:3.8125)N1190:6.8125,(Carex_tomentosa:12.300001,(Carex_melanostachya:10.2,(Carex_pseudocyperus:9.439999,(Carex_distans:8.1,(Carex_pendula:5.3,Carex_demissa:5.3)N1379:2.8)N1370:1.339999)N1326:0.76)N1313:2.100001)N1256:21.199999)N1185:40.024992)N1063:74.27501,(Chelidonium_majus:135,((((Hypericum_montanum:5,Hypericum_hirsutum:5)N2928:98.099997,Astragalus_glycyphyllos:103.099998)N2698:5.900001,((Descurainia_sophia:27.299999,((Lepidium_campestre:13.2,Coronopus_squamatus:13.2)Lepidieae:12.373334,(Camelina_sativa:6.6,Camelina_microcarpa:6.6,Camelina_alyssum:6.6)Camelina:18.973334)N4741:1.726666)Lineage_I_Franzke_et_al._2011:5,Bunias_orientalis:32.300001,((Berteroa_incana:19.3,Alyssum_alyssoides:19.300001)Alysseae:12.325,Arabis_hirsuta:31.625,(Isatis_tinctoria:29.600002,Diplotaxis_tenuifolia:29.6)Lineage_II_Franzke_et_al._2011:2.024999)Expanded_Lineage_II_Franzke_et_al._2011:0.675)Core_Brassicaceae:76.700001)Rosidae:9.199999,((((Dianthus_deltoides:4.53,Dianthus_gratianopolitanus:4.53)Dianthus:27.04,Holosteum_umbellatum:31.569999)N5370:27.755001,(Bassia_scoparia:51.730001,Atriplex_rosea:51.73)N5582:7.594999)N5289:52.175002,Adenostyles_alliariae:111.500003)N5136:6.7)N2487:16.800009)Eudicots:12.800003)N398;"

testContMap(rawBadTree)
testNodes(rawBadTree)

testContMap(rawGoodTree)
testNodes(rawGoodTree)

If I run the testContMap with rawGoodTree, it works properly. If I run the testContMap with rawBadTree, it reports error.

It seems that the root cause is the following use of matchNodes. It returns NAs for the rawBadTree: https://github.com/liamrevell/phytools/blob/75de27a0eb2e06682f7634dd3fa47098a8d408dd/R/fastAnc.R#L37

testNodes(rawGoodTree)
    tr1 tr2
[1,]  34  34
[2,]  35  35

testNodes(rawBadTree)
     tr1 tr2
[1,]  35  NA
[2,]  36  NA

How can I use contMap function with rawBadTree? Thanks for your help.

[liamrevell]

After spending some time on this it appears that the problem is with ape::multi2di. You should be able to see this by running the following example:

library(ape)
library(phytools)
bad <- "((1:73.52,((2:26.69,((3:5.72,4:5.72):17.16,5:22.88):3.81):6.81,(6:12.3,(7:10.2,(8:9.44,(9:8.1,(10:5.3,11:5.3):2.8):1.34):0.76):2.1):21.2):40.02):74.28,(12:135,((((13:5,14:5):98.1,15:103.1):5.9,((16:27.3,((17:13.2,18:13.2):12.37,(19:6.6,20:6.6,21:6.6):18.97):1.73):5,22:32.3,23:32.3,((24:19.3,25:19.3):12.32,26:31.62,(27:29.6,28:29.6):2.02):0.68):76.7):9.2,((((29:4.53,30:4.53):27.04,31:31.57):27.76,(32:51.73,33:51.73):7.59):52.18,34:111.5):6.7):16.8):12.8);"
tt<-read.tree(text=bad)
bt<-multi2di(tt)
par(mfrow=c(1,2))
plotTree(tt,fsize=0.7)
plotTree(bt,fsize=0.7,direction="leftwards")

Here's the good news though. Evidently this is already updated in ape version 5.6 -- so all you need to do is update ape from CRAN and you'll be set.