I had a look at the tree. There seems a bug in extract.clades, when the
tree was rooted before with the root() function. extract.clade() seem to
expect some ordering of nodes which seems not satisfied (and it looks as
it was written a long time ago).
I just wrote a small function (easier than to debug the old one) which
should do the same and it seems to work:
library(phangorn)
extract.clade2 <- function(phy, node){
if(!(node %in% phy$edge[,1])) stop("node number must be greater than
the number of tips")
drop.tip(phy, setdiff(1:Ntip(phy), Descendants(phy, node)[[1]]))
}
Klaus
On Jan 5, 2017 10:25 PM, "Joseph W. Brown" <josep...@umich.edu
<mailto:josep...@umich.edu>> wrote:
Hmm. Maybe something wonky with your tree? I simulated a 17-tip tree
and tried things out:
phy <- rtree(17);
rootID <- length(phy$tip.label) + 1;
counter <- 1;
for (i in rootID:max(phy$edge[,1])) {
clade <- extract.clade(phy, i);
# do something. just printing clade properties here
print(paste0("clade ", counter, " (node ", i, ") has ",
length(clade$tip.label), " tips."));
counter <- counter + 1;
}
[1] "clade 1 (node 18) has 17 tips."
[1] "clade 2 (node 19) has 11 tips."
[1] "clade 3 (node 20) has 10 tips."
[1] "clade 4 (node 21) has 4 tips."
[1] "clade 5 (node 22) has 3 tips."
[1] "clade 6 (node 23) has 2 tips."
[1] "clade 7 (node 24) has 6 tips."
[1] "clade 8 (node 25) has 5 tips."
[1] "clade 9 (node 26) has 3 tips."
[1] "clade 10 (node 27) has 2 tips."
[1] "clade 11 (node 28) has 2 tips."
[1] "clade 12 (node 29) has 6 tips."
[1] "clade 13 (node 30) has 2 tips."
[1] "clade 14 (node 31) has 4 tips."
[1] "clade 15 (node 32) has 3 tips."
[1] "clade 16 (node 33) has 2 tips."
and had no problem. Maybe post your tree here?
JWB
________________________________________
Joseph W. Brown
Post-doctoral Researcher, Smith Laboratory
University of Michigan
Department of Ecology & Evolutionary Biology
Room 2071, Kraus Natural Sciences Building
Ann Arbor MI 48109-1079
josep...@umich.edu <mailto:josep...@umich.edu>
> On 5 Jan, 2017, at 17:18, Kamila Naxerova <knaxer...@partners.org
<mailto:knaxer...@partners.org>> wrote:
>
> Dear Joseph,
>
> thanks so much. This is exactly what I need!
>
> I am running into some problems that I don’t understand though. In
my case, rootID is 18, and max(phy$edge[,1]) is 33. When I try to
execute your loop, this happens:
>
> > extract.clade(phy, 18)
>
> Phylogenetic tree with 17 tips and 16 internal nodes.
>
> Tip labels:
> X1, X8, X9, X10 ...
>
> Rooted; includes branch lengths.
>
> So far so good… but then I keep getting these errors:
>
> > extract.clade(phy, 19)
> Error in phy$edge[, 2] : incorrect number of dimensions
> > extract.clade(phy, 20)
> Error in phy$edge[, 2] : incorrect number of dimensions
>
> Not sure why extract.clade produces these errors. 19-23 don’t
work, 24-26 work, 27 produces the error again, 28 works etc.
>
> Thanks again.
>
> Kamila
>
>
>> On Jan 5, 2017, at 4:12 PM, Joseph W. Brown <josep...@umich.edu
<mailto:josep...@umich.edu> <mailto:josep...@umich.edu
<mailto:josep...@umich.edu>>> wrote:
>>
>> Not sure if I understand the problem completely, but this should
allow you to examine all of the clades (and should work if
polytomies are involved):
>>
>> # for tree phy
>> rootID <- length(phy$tip.label) + 1;
>> for (i in rootID:max(phy$edge[,1])) {
>> clade <- extract.clade(phy, i);
>> # do something
>> }
>>
>> This includes the root node (i.e. whole tree), but that can be
changed. This can be rewritten as an lapply if necessary.
>>
>> HTH.
>> JWB
>> ________________________________________
>> Joseph W. Brown
>> Post-doctoral Researcher, Smith Laboratory
>> University of Michigan
>> Department of Ecology & Evolutionary Biology
>> Room 2071, Kraus Natural Sciences Building
>> Ann Arbor MI 48109-1079
>> josep...@umich.edu <mailto:josep...@umich.edu>
<mailto:josep...@umich.edu <mailto:josep...@umich.edu>>
>>
>>
>>
>>> On 5 Jan, 2017, at 15:50, Kamila Naxerova
<knaxer...@partners.org <mailto:knaxer...@partners.org>
<mailto:knaxer...@partners.org <mailto:knaxer...@partners.org>>>
wrote:
>>>
>>> Hi all,
>>>
>>> I would like to break a phylogenetic tree into all possible
clades and then examine each one of them for certain characteristics.
>>>
>>> I am writing some code to do this from scratch, but it’s getting
pretty cumbersome quickly.
>>>
>>> I was wondering whether there are some functions out there
already that could help me with this task?
>>>
>>> Thanks so much for any help.
>>>
>>> Cheers,
>>> Kamila
>>>
>>>
>>> The information in this e-mail is intended only for the person
to whom it is
>>> addressed. If you believe this e-mail was sent to you in error
and the e-mail
>>> contains patient information, please contact the Partners
Compliance HelpLine at
>>> http://www.partners.org/complianceline
<http://www.partners.org/complianceline>
<http://www.partners.org/complianceline
<http://www.partners.org/complianceline>> . If the e-mail was sent
to you in error
>>> but does not contain patient information, please contact the
sender and properly
>>> dispose of the e-mail.
>>> _______________________________________________
>>> R-sig-phylo mailing list - R-sig-phylo@r-project.org
<mailto:R-sig-phylo@r-project.org> <mailto:R-sig-phylo@r-project.org
<mailto:R-sig-phylo@r-project.org>>
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
<https://stat.ethz.ch/mailman/listinfo/r-sig-phylo>
<https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
<https://stat.ethz.ch/mailman/listinfo/r-sig-phylo>>
>>> Searchable archive at
http://www.mail-archive.com/r-sig-phylo@r-project.org/
<http://www.mail-archive.com/r-sig-phylo@r-project.org/>
<http://www.mail-archive.com/r-sig-phylo@r-project.org/
<http://www.mail-archive.com/r-sig-phylo@r-project.org/>>
>
[[alternative HTML version deleted]]
_______________________________________________
R-sig-phylo mailing list - R-sig-phylo@r-project.org
<mailto:R-sig-phylo@r-project.org>
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
<https://stat.ethz.ch/mailman/listinfo/r-sig-phylo>
Searchable archive at
http://www.mail-archive.com/r-sig-phylo@r-project.org/
<http://www.mail-archive.com/r-sig-phylo@r-project.org/>