-
Notifications
You must be signed in to change notification settings - Fork 57
Description
Hi Liam! I hope you are enjoying the holidays.
I wrote a short script to visualize GC content on a tree. I need to use external sequence reconstructions, and I have:
- a newick file with names like
A30on the internal nodes - a FASTA file with sequences for the internal nodes.
It looks like contMap( ) crashes if anc.states contains node names:
> contMap(tree, gc_values[tree$tip.label],method="user",anc.states=gc_values[tree$node.label])
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
I was able to make it work by doing the following:
anc.states = gc_values[tree$node.label]
names(anc.states) <- NULL
contMap(tree, gc_values[tree$tip.label],method="user",anc.states=anc.states)
It seems like the code is expecting that any ancestral sequences have names that are just integers represented as strings, even if node.label has other names:
else {
if (is.null(names(anc.states)))
names(anc.states) <- 1:tree$Nnode + Ntip(tree)
a <- anc.states[as.character(1:tree$Nnode + Ntip(tree))]
}
Here, if names(anc.states) is not null, but also isn't equal to the integer-as-string names, then doing a <- anc.states[as.character(1:tree$Nnode + Ntip(tree))] creates a vector a that is filled with NA values.
Also, I found it a bit confusing that I had to split the values into different vectors for the tip nodes and the ancestral nodes.