xavierdidelot / TransPhylo

Reconstruction of transmission trees using genomic data
http://xavierdidelot.github.io/TransPhylo/
GNU General Public License v2.0
60 stars 22 forks source link

Errors when using infer_multittree_share_param #32

Closed lingxuan85511 closed 10 months ago

lingxuan85511 commented 10 months ago

Dear author,

I am usig infer_multittree_share_param (TransPhylo v1.4.10) to infer the parameters of transmission chain with my data. I found this function performed well when the dataset is not big. However, error occurred when I used a big dataset(>3000 tips). I have no idea how to solve this problem. Do you have any suggestions? The script and error information was as follow:

library(ape) library(TransPhylo) library(coda) packageVersion("TransPhylo") [1] ‘1.4.10’ tree1<-multi2di(read.nexus("timetree_21.nexus")) tree1$edge.length[tree1$edge.length<=0]<-0.1/365 tree1$edge.length<-tree1$edge.length*365 ptree1<-ptreeFromPhylo(tree1,dateLastSample=60)

tree2<-multi2di(read.nexus("timetree_36.nexus")) tree2$edge.length[tree2$edge.length<=0]<-0.1/365 tree2$edge.length<-tree2$edge.length*365 ptree2<-ptreeFromPhylo(tree2,dateLastSample=60)

ptrees<-list(ptree1,ptree2) w.shape<-2.892313 w.scale<-2.938824 dateT<-61 record<-infer_multittree_share_param(ptrees,mcmcIterations=1000000,thinning=100,w.shape=w.shape,w.scale=w.scale,updateNeg=T,updateOff.p=T,dateT=dateT,share=c("off.r","off.p","pi"),delta_t=1)

|
| | 0%Error in if (log(runif(1)) < sum(pTTree2) - sum(pTTree)) { : missing value where TRUE/FALSE needed Calls: infer_multittree_share_param ... with -> with.default -> eval -> eval -> one_update_share Execution halted

xavierdidelot commented 10 months ago

Does it work if you use only ptree1 or ptree2?

lingxuan85511 commented 10 months ago

Yese, when using inferTTree function with each tree, the program is running correctly.

xavierdidelot commented 10 months ago

Ok, in that case the problem is specific to the infer_multitree_share_param that was written by Yuanwei Xu. I don't think he monitors the issues on this repository though, so you should probably try to contact him directly using his email address provided in the DESCRIPTION file.

lingxuan85511 commented 10 months ago

Thank you for your kindly help. I will directly contact him via email.