problem with Chainsaw2

[stiegosaur]

Hi, I've been using Chainsaw2 outside of Biogeobears, and I've noticed that sometimes when I make particularly big cuts (~100 million years or more on a tree with the root at 250 MA) I end up with a tree that has a lot of negative branches (see attached file). Sometimes I can get around this by making several cuts in a row, but even that is not 100% reliable which is especially troublesome inside of a loop. If anyone else has encountered this or can think of reasons why this might be happening it would be helpful.

[Nick Matzke]

Eeek!  If you email the code/example file, I can take a look.  Cheers, Nick

On Thu, Sep 24, 2015 at 1:57 PM, stiegosaur <stieg...@gmail.com> wrote:

Hi, I've been using Chainsaw2 outside of Biogeobears, and I've noticed that sometimes when I make particularly big cuts (~100 million years or more on a tree with the root at 250 MA) I end up with a tree that has a lot of negative branches (see attached file). Sometimes I can get around this by making several cuts in a row, but even that is not 100% reliable which is especially troublesome inside of a loop. If anyone else has encountered this or can think of reasons why this might be happening it would be helpful.

--
You received this message because you are subscribed to the Google Groups "BioGeoBEARS" group.
To unsubscribe from this group and stop receiving emails from it, send an email to biogeobears...@googlegroups.com.
To post to this group, send email to bioge...@googlegroups.com.
Visit this group at http://groups.google.com/group/biogeobears.
To view this discussion on the web visit https://groups.google.com/d/msgid/biogeobears/08a50774-8862-4bce-8169-f1d53cc03188%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[stiegosaur]

Thanks. I took mrbayes trees and ran them through Nylander's burntree perl script to scale the branch lengths.  The simplest code (below) yielded the attached jpeg using the attached "a.tree". The trees come out of the perl script unladderized, and I've tried ladderizing on the right and left sides, saving the ladderized trees to a new file and reading them back in, etc.

tree<-read.tree("a.tree")

x<-chainsaw2(tree,timepoint=150)

x<-x$tree_to_chainsaw

plot(x,cex=.01)

[Nick Matzke]

OK, I can replicate the issue, but I'm not sure I can fix it right away -- it's probably something about the fossil branches, which require special handling when modifying the R tree object when chainsawing (imagine temporarily extending the fossil tips to present, then chopping down the fake extended branches, then chopping back to the original branch heights, minus the bit chainsawed off the top of the tree. etc.).  It may be that this strategy backfires if big-enough cuts are made (e.g., whole clades of fossil-only taxa are being cut).

Anyhoo -- I found that it works if you just run a loop of small progressive chainsaw cuts.  You could then save just some of these, or whatever.  So that's a workaround for now...

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

#######################################################

# chainsaw2 issue

#######################################################

# https://groups.google.com/forum/?utm_medium=email&utm_source=footer#!msg/biogeobears/Jy9uYckOL7s/E3lsEJHgAwAJ

library(BioGeoBEARS)

#source("/drives/Dropbox/_njm/__packages/BioGeoBEARS_setup/BioGeoBEARS_stratified_v1.R")

source(http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_stratified_v1.R)

# Change to your working directory

wd = "/drives/Dropbox/_njm/__packages/BioGeoBEARS_setup/z_help/chainsaw2/"

setwd(wd)

#######################################################

# I'm not sure why this problem happens -- probably 

# something involving the fossils.  Fossils require

# weird contortions in the R phylo tree structure,

# when chainsaw-ing the tree.

#

# (E.g., perhaps the issue is if the time-slice is 

#  below a whole clade of fossil-only tips)

#######################################################

tree<-read.tree("a.tree")

plot(tree)

axisPhylo()

x<-chainsaw2(tree,timepoint=150)

x<-x$tree_to_chainsaw

plot(x,cex=.01)

#######################################################

# Anyway, it works if you do a series of small time slices:

#######################################################

pdffn = "chainsaw_results.pdf"

pdf(file=pdffn, width=8.5, height=22)

timecuts = seq(0.5, 249, 0.5)

tmptree = tree

new_tree_top_absolute = 0

new_tree_top_relative = 0

cat("\nChainsawing at: ")

for (i in 1:length(timecuts))

{

cat(timecuts[i], " ", sep="")

new_tree_top_relative = timecuts[i] - new_tree_top_absolute

x = chainsaw2(tmptree, timepoint=new_tree_top_relative)

tmptree = x$tree_to_chainsaw

plot(tmptree,cex=.01)

axisPhylo2(side=1, roundlabels=TRUE, minage=timecuts[i])

new_tree_top_absolute = timecuts[i]

}

dev.off()

cmdstr = paste0("open ", pdffn)

system(cmdstr)

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

[Nick Matzke]

Ah, so it dies when the cut goes beyond the bottom of the (only) branch extending to time=0:

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

tree<-read.tree("a.tree")

plot(tree, show.tip.label=FALSE)

axisPhylo()

trtable = prt(tree, printflag=FALSE)

tree_height = max(trtable$node_ht)

abline(v=tree_height-127.5, lty="dashed", col="grey")

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

Not sure how to fix it, but the progressive-trimming strategy works...

Cheers!

Nick

--

You received this message because you are subscribed to the Google Groups "BioGeoBEARS" group.
To unsubscribe from this group and stop receiving emails from it, send an email to biogeobears...@googlegroups.com.
To post to this group, send email to bioge...@googlegroups.com.
Visit this group at http://groups.google.com/group/biogeobears.

To view this discussion on the web visit https://groups.google.com/d/msgid/biogeobears/90021b24-98ff-4c0a-b3e7-25059e3b09c0%40googlegroups.com.

[stiegosaur]

Ok, thanks a bunch.

On Wednesday, September 23, 2015 at 11:57:16 PM UTC-4, stiegosaur wrote: