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> 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



> On 5 Jan, 2017, at 17:18, Kamila Naxerova <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>> 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>
>>
>>
>>
>>> On 5 Jan, 2017, at 15:50, Kamila Naxerova <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> . 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>
>>> 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-ph...@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
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at 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
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

Reply via email to