BEAST result in R

[Juliet Lindo]

Hi,

I trying to analyze my divergence time analysis result from BEAST in R, and I encounter this error, 

Error in attr(data, "mcpar") <- c(start, end, thin) : 

  attempt to set an attribute on NULL

I hope someone can help me...

Thank you very much,

Juliet

[Tim Vaughan]

Hi Juliet, what R packages are you using?  Which BEAST data files are you trying to load into R?

Tim

--
You received this message because you are subscribed to the Google Groups "beast-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to beast-users...@googlegroups.com.
To post to this group, send email to beast...@googlegroups.com.
Visit this group at https://groups.google.com/group/beast-users.
For more options, visit https://groups.google.com/d/optout.

[Santiago Sánchez]

Hi Juliet,

What packages are you using to load your BEAST data? That I know of, there is phyloch and OutBreakTools, which work pretty well. phyloch uses the function read.beast to generate phylo-object-associated vectors for all the metadata in the tree (e.g. posterior, 95%HPD, rates, etc). Here the data is associated to the node numbering and order. OutBreakTools uses read.annotated.nexus to do fairly the same, but the values are associated with edges (numbering and order), not nodes.

Hope this helps,

Santiago

--

You received this message because you are subscribed to the Google Groups "beast-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to beast-users...@googlegroups.com.
To post to this group, send email to beast...@googlegroups.com.
Visit this group at https://groups.google.com/group/beast-users.
For more options, visit https://groups.google.com/d/optout.

--

==========================

Santiago Sanchez-Ramirez, PhD

Postdoctoral Associate

Ecology and Evolutionary Biology

University of Toronto

==========================

[Juliet Lindo]

Hi Tim,

Im using these packages, 

> library(ape)

> library(maps)

> library("phytools")

> library("colorspace")

> library("XML")

> library("phyloch")

> library("geoscale")

> library("strap")

> library("coda")

and the file I'm analyzing is the .tre file, my summary tree from Tree Annotator.

~ Juliet 

[Tim Vaughan]

Thanks Juliet, which R command are you actually using to to load in the data?

--

[Juliet Lindo]

Hi Tim and Santiago,

Im using the script as a guide from the sample file I got form the net, 

Here's the script

library("phytools")

library("phyloch")

library("strap")

library("coda")

t <- read.beast("bearsDivtime_FBD.summary.tre")

t$root.time <- t$height[1]

log_data1 <- read.table("bearsDivtime_FBD.1.log",header=T)

log_data2 <- read.table("bearsDivtime_FBD.2.log",header=T)

nMCMC <- length(log_data1$originFBD)-1

burnin <- 0.2

id1 <- as.integer(nMCMC*burnin)+1

id2 <- nMCMC+1

comb_origin_data <- c(log_data1$originFBD[id1:id2],log_data2$originFBD[id1:id2])

origin_mcmc <-as.mcmc(comb_origin_data)

origin_time <- mean(origin_mcmc)

stem_length <- origin_time - t$root.time

origin_HPD <- HPDinterval(origin_mcmc)

num_taxa <- length(t$tip.label)

display_all_node_bars <- FALSE

names_list <-vector()

for (name in t$tip){

  v <- strsplit(name,"_")[[1]]

  if(display_all_node_bars){

  names_list = c(names_list,name)

  }

  else if(v[length(v)]=="0"){

  names_list = c(names_list,name)

  }

}

nids <- vector()

pos <- 1

len_nl <- length(names_list)

for(n in names_list){

  for(nn in names_list[pos:len_nl]){

    if(n != nn){

      m <- getMRCA(t,c(n,nn))

      if(m %in% nids == FALSE){

        nids <- c(nids,m)

      }

    }

  }

  pos<-pos+1

}

for(tp in 1:length(t$tip.label)){

  v <- strsplit(t$tip.label,"_")[[tp]]

  new_l <- v[1]

  if(length(v) > 2){

    new_l <- paste(v[1],v[2],sep="_")

  }

  if(tail(v,1) != "0"){

  new_l <- paste(new_l, "X",sep=" ")

  }

  t$tip.label[tp] <- new_l

}

root_max <- t$"height_95%_HPD_MAX"[1]

x_max <- origin_HPD[2] * 0.1 + origin_HPD[2]

pdf("geoscaled_bears.pdf", width=10, height=7)

geoscalePhylo(tree=ladderize(t,right=FALSE), boxes="Age", cex.tip=1.2,cex.age=1,

cex.ts=1,label.offset=0,x.lim=c(-15,x_max),lwd=1.5)

lastPP<-get("last_plot.phylo",envir=.PlotPhyloEnv)

origin_xx <- c(lastPP$xx[num_taxa+1],-stem_length)

lines(origin_xx,c(lastPP$yy[num_taxa+1],lastPP$yy[num_taxa+1]),lwd=1.5,col="black")

bar_xx_o <- c(t$root.time-origin_HPD[1], t$root.time-origin_HPD[2])

lines(bar_xx_o,c(lastPP$yy[num_taxa+1],lastPP$yy[num_taxa+1]),col=rgb(1,0,0,alpha=0.3),lwd=8)

points(-stem_length,lastPP$yy[num_taxa+1],pch=15,col="red",cex=1.5)

for(nv in nids){

  bar_xx_a <- c(lastPP$xx[nv]+t$height[nv-num_taxa]-t$"height_95%_HPD_MIN"[nv-num_taxa], lastPP$xx[nv]-(t$"height_95%_HPD_MAX"[nv-num_taxa]-t$height[nv-num_taxa]))

  lines(bar_xx_a,c(lastPP$yy[nv],lastPP$yy[nv]),col=rgb(0,0,1,alpha=0.3),lwd=8)

}

t$node.label<-t$posterior

p <- character(length(t$node.label))

p[t$node.label >= 0.95] <- "black"

p[t$node.label < 0.95 & t$node.label >= 0.75] <- "gray"

p[t$node.label < 0.75] <- "white"

nodelabels(pch=21, cex=1.5, bg=p)

dev.off()

On Tuesday, March 28, 2017 at 8:34:14 PM UTC-4, Juliet Lindo wrote:

[Juliet Lindo]

Hi Tim,

But when I get to this part

"origin_mcmc <-as.mcmc(comb_origin_data)", I get the error.

I edit some of the scripts for my data. Prior the to command above, all commands work,

~ Juliet

On Tuesday, March 28, 2017 at 8:34:14 PM UTC-4, Juliet Lindo wrote:

[Tim Vaughan]

Hi Juliet, can you confirm that the data frames "log_data1" and "log_data2" contain what they are meant to after loading in the log files? Specifically, do they contain columns named originFBD and are do those columns contain finite (non NaN or NA) values?

--

[Juliet Lindo]

Hello  Tim,

I've checked the sample .log file, and one of them contain this

 <parameter id="originFBD.t:bearsTree" lower="0.0" name="origin">45.5</

Any of my .log files dont have information. Is it because when running BEAST in the sample file, Fossilised Birth Death Model was used, while for my analysis, I used Yule Model?

Thank you so much for the help.

~ Juliet

[Juliet Lindo]

I was trying to view my tree in R, I know I can do that in Figtree but the features and info in it is limited.  In R, I can have a better figure because i can add the node bars, ages and geologic time scale...

~ Juliet

[Tim Vaughan]

Hi Julia, sorry - when you said that you had changed the script to suit your data, I thought you'd also changed lines like "comb_origin_data <- c(log_data1$originFBD[id1:id2],log_data2$originFBD[id1:id2])" which are absolutely dependent on the data file.  At very least you'll need to change these references to correspond to parameters in the model you're using.  I don't know whether the rest of this particular script is translatable in this way, but there are likely other changes you'll have to make as well.

Hope this helps a little bit,

Tim

[Juliet Lindo]

Hi Tim,

Its a huge help. At least now I know where the error originated from and now I can figure out what problem should I address. 

Thank you very much,

Juliet

[Santiago Sánchez]

Juliet,

Maybe just to add a bit. To plot your tree and visualize summarized data on it, you only need to load your MCC tree, then plot and add info to it.

For example:

library(phyloch)

t <- read.beast("bearsDivtime_FBD.summary.tre")

plot(ladderize(t))

HPDbars(t)

add.geoscale(t)

This will give you just a basic plot. Of course you can tweak it to look nicer.

Have a look at the examples of each one of these function to give you ideas.

I'm no sure what the final end-purpose of the script is, but, again, if you only wish to visualize your tree, you can go for something simpler.

Cheers,

Santiago

--

[Juliet Lindo]

Hi Santiago,

Thank you for the info. Its a big help for me. Appreciate it.

~ Juliet

[Juliet Lindo]

Hi Tim and Santiago,

I was able to work it out. I have my tree now from R. Thank you so much for the help.

~ Juliet

[Deborah Veranea Espinosa Martinez]

Hi Juliet

My name is Deborah and I want to use this script, is it possible? and Where I can check the command description?

Thaks 

Deborah

[HS]

Dear Juliet,

I would suggest to use ggtree package https://yulab-smu.github.io/treedata-book/index.html

It is a good one, has a lot of smart features. I believe with this package you can diminish the size of your script profoundly. Also, you can have many nice features.

Best,

Hovhannes

[Alexander Sáchez Rojas]

Hi Juliet, I would be very grateful if you share the solution of this error. I've run into this same problem and been flipping the code for a long time, but I can't find a solution.

I hope you can help me

Thank you very much,
Alex