Re: [R-sig-phylo] chronos ape substitution rate strict clock

[Klaus Schliep]

Dear Vincenzo,
For ML estimates the edge length are the expected number of substitutions per 
site, which depends on the product of rate and time. So you need either a rate 
estimate or calibration points to estimate the time. If you divide the edge 
length by the rate the edges should be proportional to time.
With pml_bb you can estimate ultrametric or tip dated trees with a strict 
molecular clock model.
Kind regards,
Klaus




Von meinem iPhone gesendet

> Am 04.03.2024 um 21:06 schrieb Vincenzo Ellis :
> 
> Dear Liam,
> 
> Thank you for your answer. Yes, I have a hypothesized average clock rate
> and no explicit calibration points.
> 
> If I use ape::chronos() with default values it gives me a depth of 1 at the
> root. So I suppose to rescale the branch lengths I just need to multiply
> all of the branch lengths by a value (i.e., tree$edge.length * value), is
> that right?
> 
> Thanks again,
> 
> Vincenzo
> 
>> On Mon, Mar 4, 2024 at 11:43 AM Liam J. Revell  wrote:
>> 
>> Dear Vincenzo.
>> 
>> If I understand your problem problem, you do not have any explicit
>> calibration points -- but you have a hypothesized average clock rate?
>> 
>> If so, then you can obtain an ultrametric tree from* ape::chronos* for
>> any value of the smoothing parameter (*lambda*) and then simply re-scale
>> it to have the desired total depth (based on your hypothesized clock rate).
>> 
>> To choose a "correct" value of *lambda* one can use cross-validation as
>> described in Sanderson (2002; doi:10.1093/oxfordjournals.molbev.a003974).
>> 
>> Others should feel welcome to weigh in if this is not right.
>> 
>> All the best, Liam
>> Liam J. Revell
>> Professor of Biology, University of Massachusetts Boston
>> Web: http://faculty.umb.edu/liam.revell/
>> Book: Phylogenetic Comparative Methods in R
>> 
>> (*Princeton University Press*, 2022)
>> 
>> 
>> On 3/4/2024 11:09 AM, Vincenzo Ellis wrote:
>> 
>> [You don't often get email from vael...@udel.edu. Learn why this is 
>> important at https://aka.ms/LearnAboutSenderIdentification ]
>> 
>> CAUTION: EXTERNAL SENDER
>> 
>> Dear Emmanuel,
>> 
>> Thank you very much for your response. I cannot see how to provide the
>> substitution rate to the phangorn::pml_bb() function, but I was looking at
>> the ape::node.dating() function and it appears that I could provide the
>> substitution rate to the "mu" argument and then set the "node.dates"
>> argument to NA or zero for all tips (I'm not sure if NA or zero would be
>> preferable to force the tips to all be from a single time point). Do you
>> think that would work? I'm not sure how to make ape::node.dating() accept a
>> substitution rate rather than try to estimate one. Maybe an option could be
>> added to allow mu to equal a user-specified number rather than the output
>> of ape::estimate.mu()?
>> 
>> Another option might be to calculate an estimated age for every node
>> connecting sister taxa in the tree by converting the genetic distances
>> between sister pairs to divergence times using the substitution rate and
>> then use those as priors in ape::chronos(). I suppose I could also apply
>> that logic to date all of the nodes by using the mean pairwise distances
>> between taxa on either side of a node and converting that to divergence
>> times (although the R code for such a calculation would likely take me a
>> while to figure out). Would that be another option?
>> 
>> Thanks again,
>> 
>> Vincenzo
>> 
>> On Sat, Mar 2, 2024 at 5:44 AM Emmanuel Paradis  
>> 
>> wrote:
>> 
>> 
>> Hi Vincenzo,
>> 
>> There's no direct way to do this with ape::chronos(). You may have a look
>> at the function phangorn::pml_bb() but I'm not sure it can estimate the
>> dates if the rate is provided in a model object given as main argument(?)
>> 
>> That said, I expect that estimating so many dates to be very challenging
>> (unless you have a lot of known dates for calibration). This implies that
>> you are certainly right to look for an approach where you don't need to
>> estimate the rates.
>> 
>> Best,
>> 
>> Emmanuel
>> 
>> - Le 27 Fév 24, à 22:40, Vincenzo Ellis vael...@udel.edu a écrit :
>> 
>> 
>> Dear R-sig-phylo members,
>> 
>> I've made a maximum likelihood tree in Raxml for several thousand taxa
>> using a single gene that has an estimated substitution rate of 0.006
>> substitutions/nucleotide/My. Is there a way to use chronos in ape to
>> 
>> apply
>> 
>> that substitution rate as a fixed clock rate and generate an ultrametric
>> time-scaled version of the tree?
>> 
>> Thank you,
>> 
>> Vincenzo
>> 
>>  [[alternative HTML version deleted]]
>> 
>> ___
>> R-sig-phylo mailing list - 
>> 

Re: [R-sig-phylo] Parallelization in ape::dist.topo

[Klaus Schliep]

Dear Vojtěch,
nice work. Just a few random comments:
Parallelization is often not straightforward as it depends on the hardware
and the operating system. My preference is using the future package for
parallelization as it does some nice abstraction for the different R
packages, so you can try different things and no need to write different
versions of the code. chatGPT seems to have missed this.
Defaulting to "cores = detectCores()" is not a good idea. The CRAN
Repository Policy (https://cran.r-project.org/web/packages/policies.html)
states: "If running a package uses multiple threads/cores it must never use
more than two simultaneously: the check farm is a shared resource and will
typically be running many checks simultaneously." In short such a default
can cause serious problems on a cluster and I got properly ripleyed for
this years ago, but that's another story. So the default should be
something like min(2L, detectCores()). So with great power comes great
responsibility and users should read the man pages.
Kind regards,
Klaus


On Mon, Mar 6, 2023 at 2:43 PM Vojtěch Zeisek  wrote:

> Hello dear colleagues,
> I use often ape::dist.topo (see here dist.topo.r), which is doing the
> calculations sequentially, which is very slow for large data sets.


What is large? Large number of trees or large trees?

I'm sorry,
> I haven't found any relevant Git repository or so, so I hope Emmanuel
> won't
> mind if I discuss it here.
> I discussed various options with ChatGPT and dist.topo.par1.r is the
> simplest
> solution, basically using mc.lapply instead of 2 for loops. Good study
> material for how to do it in general. Little enhancements are in
> dist.topo.par2.r, which should be slightly better in case some pair of
> comparisons would return NA or so, but from my tests there doesn't seem to
> be
> any difference.
>

I think there is room for improvement with preprocessing trees. If you have
e.g. bootstrap trees with short edges you might want to use di2multi()
first to avoid spurious differences. Filtering out duplicated trees could
also speed up things, if there are many. This will depend on the trees you
want to compare and the methods.

And finally there is dist.topo.par3.r which doesn't load parallel (and uses
> plain lapply) for cores==1, while parallel and doParallel for multiple
> cores.
> It also contains some checks and error handling. From my testing it works
> well. I'm not sure if tryCatch is really needed there.


I am not sure if it is necessary, but you don't want to have it in the
inner loop ;)

In any case,
> improvements welcomed. :-)
> So, what do You think? Is this usable improvement of ape::dist.topo?
> Sincerely,
> V.
>
> --
> Vojtěch Zeisek
> https://trapa.cz/en/
>
> Department of Botany, Faculty of Science
> Charles University, Prague, Czech Republic
> https://www.natur.cuni.cz/biology/botany/
> https://lab-allience.natur.cuni.cz/
>
> Institute of Botany, Czech Academy of Sciences
> Průhonice, Czech Republic
> https://www.ibot.cas.cz/en/
> Computing cluster
> https://sorbus.ibot.cas.cz/en/start
> ___
> 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/
>


-- 
Klaus Schliep

Senior Scientist
Institute of Molecular Biotechnology
TU Graz
https://www.imbt.tugraz.at

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

Re: [R-sig-phylo] ancestral state reconstruction with a complex discrete, trait

[Klaus Schliep]

Dear list,
sorry for being late to the discussion. Thanks to the generous help of Iris
Bardel-Kahr we are improving the vignettes for phangorn adding some more
examples for morphological / complex / special traits. Please have a look
at the vignettes on the github page https://klausvigo.github.io/phangorn/
in the articles section. The titles might change and sections still get
moved around at the moment. So what you are interested in is distributed
over several vignettes. And suggestions for improvement are welcome.
In optim.pml() / pml_bb() you need to change the default value to
optEdge=FALSE in case you only want to estimate the parameters for the rate
matrix. Only time reversible models are possible at the moment.
A parsimony analysis at the start might be useful. If the parsimony score
divided by the number of sites is high, if you expect many state changes
the ancestral reconstruction. At the root you will see a mixture of all
states, which is a sign of uncertainty. This can also be an artefact if you
observe transitions on short edges (so-called cherries) which will bias
transition rates. However this you can catch if the number of expected
substitutions is considerably higher than the parsimony score / site.
Kind regards,
Klaus

On Mon, Nov 14, 2022 at 8:03 AM Jonas Eberle  wrote:

> Dear Chris,
>
> apart from the phangorn citation
> (https://cran.r-project.org/web/packages/phangorn/citation.html) you
> might be interested in this one:
> Sabatinelli, G., Eberle, J., Fabrizi, S. & Ahrens, D. (2020) A molecular
> phylogeny of Glaphyridae (Coleoptera: Scarabaeoidea): evolution of
> pollination and association with “Poppy guild” flowers. /Systematic
> Entomology/ n/a. https://doi.org/10./syen.12429
>
> Good luck with your study!
>
> Best,
> Jonas
>
> Am 12.11.2022 um 12:00 schrieb r-sig-phylo-requ...@r-project.org:
> > Send R-sig-phylo mailing list submissions to
> >   r-sig-phylo@r-project.org
> >
> > To subscribe or unsubscribe via the World Wide Web, visit
> >   https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
> > or, via email, send a message with subject or body 'help' to
> >   r-sig-phylo-requ...@r-project.org
> >
> > You can reach the person managing the list at
> >   r-sig-phylo-ow...@r-project.org
> >
> > When replying, please edit your Subject line so it is more specific
> > than "Re: Contents of R-sig-phylo digest..."
> >
> > Today's Topics:
> >
> > 1. Re: R-sig-phylo Digest, Vol 178, Issue 5 (Krzysztof Kozak)
> >
> > ___
> > R-sig-phylo mailing list
> > R-sig-phylo@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
>
> [[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/
>


-- 
Klaus Schliep

Senior Scientist
Institute of Environmental Biotechnology
TU Graz
https://ubt.tugraz.at <https://icbt.tugraz.at>

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

Re: [R-sig-phylo] dist.nodes

[Klaus Schliep]

Hi Nick,

it might be useful to plot the node labels on the tree:

plot(tree, label.offset = .25)
tiplabels()
nodelabels()

Regards,
Klaus




On Tue, Sep 7, 2021 at 4:56 AM Emmanuel Paradis 
wrote:

> Hi,
>
> A tree has n terminal nodes (aka tips) and m internal nodes (aka nodes
> simply). In the edge matrix, the tips are numbered 1:n and the nodes are
> numbered (n+1):(n+m) (same than n+1:m). n and m can be found with:
>
> n <- Ntip(tree) # or length(tree$tip.label)
> m <- Nnode(tree) # or tree$Nnode
>
> Details can be found in:
>
> http://ape-package.ird.fr/misc/FormatTreeR.pdf
>
> Best,
>
> Emmanuel
>
> - Le 5 Sep 21, à 21:07, Nick Youngblut nyoun...@gmail.com a écrit :
> > For `ape::dist.nodes()`, how can one match the output matrix
> rows/columns with
> > the node IDs in the tree (eg., the tip.labels)? I cannot just use
> > `ape::cophenetic()` in my particular situation. The docs for
> > `ape::dist.nodes()` state:
> >
> > ```
> > … in the case of dist.nodes, the numbers of the tips and the nodes (as
> given by
> > the element edge).
> > ```
> >
> > … but tree$edge doesn’t provide any direct info about which nodes are
> internal
> > vs external and how the external tip labels map to the tree$edge matrix.
> >
> > Thanks,
> > Nick
> > ___
> > 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/
>
> ___
> 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/
>


-- 
Klaus Schliep

Senior Scientist
Institute of Computational Biotechnology
TU Graz
https://icbt.tugraz.at
https://www.phangorn.org/ <http://www.phangorn.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/

Re: [R-sig-phylo] Using bind.tip and lapply on a multiPhylo object

[Klaus Schliep]

Dear Russell,
there is also a function getMRCA() in ape, which is more appropriate to use
for your problem, as it only computes the MRCA for a set of tips. This
should be much more efficient than using mrca().
Kind regards,
Klaus



On Thu, Jun 3, 2021 at 1:05 AM Russell Engelman 
wrote:

> Dear R-sig-phylo,
>
> I have been working with a mammalian phylogeny I recently downloaded from
> VertLife (http://vertlife.org/phylosubsets/). Unfortunately, the phylogeny
> is missing a large number of species, so I am trying to manually add these
> taxa to the phylogeny. I have a series of 100 trees that I am using to do
> things such as test for phylogenetic signal. I know how to use bind.tip to
> add new taxa to a single tree, but I am having more trouble with a
> multiPhylo object. I am primarily adding these taxa by placing them as
> sister to their nearest included relative (since most of them are elevated
> former subspecies), but the issue here is that in the 100 trees in the
> multiPhylo object the node representing the taxon to bind these taxa to is
> not the same across all trees due to shifting topologies.
>
> This is an example of the code I have been using, in which "tree" is the
> tree object. This works for a single 'phylo' tree but not 'multiphylo'.
>
> ```
> newtree<-lapply(tree,bind.tip,tip.label="Cercopithecus_albogularis",
> position=0.59,edge.length = 0.59,
>
> where=mrca(tree)["Cercopithecus_mitis","Cercopithecus_mitis"])
> ```
>
> Now, this code will not work, but I know exactly why: 'tree' is a
> multiPhylo object and so the 'where' argument cannot find the node for the
> terminal taxon. However, the issue is how can I tell R to repeat this
> 'where' argument for each of the 100 trees, since the node in question is
> not identical across these trees? Is there an easier way to do this than
> using the 'mrca' call for each terminal taxon? I've noticed adding a 'mrca'
> argument also increases computation time and if I am reinventing the wheel
> it would be nice to know if I am overthinking things.
>
> Sincerely,
> Russell
>
> [[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/
>


-- 
Klaus Schliep

Senior Scientist
Institute of Computational Biotechnology
TU Graz
https://icbt.tugraz.at
https://www.phangorn.org/ <http://www.phangorn.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/

Re: [R-sig-phylo] phangorn

[Klaus Schliep]

Dear all,
I just submitted a new version of phangorn to CRAN, the linux version is
already up there.
I hope this version survives the scrutiny of Prof. Ripley.
Kind regards,
Klaus

On Tue, Apr 6, 2021 at 12:23 PM Krzysztof Bartoszek 
wrote:

> Dear R-sig-phylo,
> I just read the worrisome email from CRAN that phagorn is scheduled for
> archival. Does anyone know if there is a plan concerning phangorn?
> Archiving all the listed packages would mean that a lot of software for the
> PCM community will be difficult to access.
>
> Best wishes
> Krzysztof Bartoszek
>
> Sent with ProtonMail Secure Email.
> [[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/
>


-- 
Klaus Schliep

Senior Scientist
Institute of Computational Biotechnology
TU Graz
https://icbt.tugraz.at
https://www.phangorn.org/ <http://www.phangorn.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/

Re: [R-sig-phylo] extracting sister taxon pairs from phylo object in APE R package

[Klaus Schliep]

l, is there some straightforward way for
> me
> >> to extract a list of sister taxon pairs (whether the sister taxa are two
> >> edges or an edge and a node)?
> >>
> >> I wrote to E. Paradis (the author of ape), and he said that he didn't
> know
> >> of a direct way to do this, but it may be possible through some
> >> manipulation of the mrca functions.
> >>
> >> Does anybody on this board have any suggestions for how to go about
> this?
> >>
> >> Thank you, Max Shpak
> >>
> >> --
> >> ===
> >> Max Shpak, Ph.D.
> >> Center for Systems and Synthetic Biology
> >> University of Texas at Austin
> >> Austin, TX 78712
> >>
> >>[[alternative HTML version deleted]]
> >>
> >> ___
> >> R-sig-phylo mailing list - R-sig-phylo@r-project.org
> >>
> https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-phylodata=02%7C01%7Cliam.revell%40umb.edu%7C90a597b4b91d49de99c908d79571b556%7Cb97188711ee94425953c1ace1373eb38%7C0%7C1%7C637142187929623185sdata=xiGcF2HLz%2Ffv6B1h0tKuQFmwLOfaneTjrXDgFNv4glA%3Dreserved=0
> >> Searchable archive at
> https://nam01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.mail-archive.com%2Fr-sig-phylo%40r-project.org%2Fdata=02%7C01%7Cliam.revell%40umb.edu%7C90a597b4b91d49de99c908d79571b556%7Cb97188711ee94425953c1ace1373eb38%7C0%7C1%7C637142187929623185sdata=XdPCQsFlTHWdRx8G6OgTczvvLEInBWVoNflalBdpkFY%3Dreserved=0
> >
> >
> >  [[alternative HTML version deleted]]
> >
> > ___
> > R-sig-phylo mailing list - R-sig-phylo@r-project.org
> >
> https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-phylodata=02%7C01%7Cliam.revell%40umb.edu%7C90a597b4b91d49de99c908d79571b556%7Cb97188711ee94425953c1ace1373eb38%7C0%7C1%7C637142187929623185sdata=xiGcF2HLz%2Ffv6B1h0tKuQFmwLOfaneTjrXDgFNv4glA%3Dreserved=0
> > Searchable archive at
> https://nam01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.mail-archive.com%2Fr-sig-phylo%40r-project.org%2Fdata=02%7C01%7Cliam.revell%40umb.edu%7C90a597b4b91d49de99c908d79571b556%7Cb97188711ee94425953c1ace1373eb38%7C0%7C1%7C637142187929623185sdata=XdPCQsFlTHWdRx8G6OgTczvvLEInBWVoNflalBdpkFY%3Dreserved=0
> >
> ___
> 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/
>


-- 
Klaus Schliep

Senior Scientist
Institute of Computational Biotechnology
TU Graz
http://www.phangorn.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/

Re: [R-sig-phylo] Testing tree disagreement: SH-test disappeared in ape-package

[Klaus Schliep]

Hi Sebastian,
there is an implementation of the SH-test and SOWH-test in phangorn. Send
me an mail off list if you have problems using them.
Regards,
Klaus

On Jun 8, 2018 8:50 AM, "Sebastian Y. Müller"  wrote:

Dear all,

I'd was wondering how to best perform a test to determine test
disagreement (congruence test) in R?
Additional information can be found in Planet et al (DOI:
10.1016/j.jbi.2005.08.008).

A common test seems the Shimodaira-Hasegawa Test (SH-test) which used to
be implemented in APE:
https://www.rdocumentation.org/packages/ape/versions/2.3-3/topics/sh.test

However, it seem to have disappeared in recent versions (e.g. ape 5.0).
Did this happen for a reason? Or maybe it has been replaced by a more
appropriate test?

Failing this, is anyone aware of other software to perform SH-, KH- or
SOWH-tests?

Thanks for any help!


-- 
Dr. Sebastian Müller
Department of Plant Sciences
University of Cambridge
Downing Street, Cambridge, CB2 3EA, UK
http://www.plantsci.cam.ac.uk/directory/mueller-sebastian
Phone: 0044 1223 748976

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

Re: [R-sig-phylo] importing trees from birdtree.org into R

[Klaus Schliep]

Dear Agustín,

the nexus file on dropbox does not contain any trees.
If you open it in an editor you see this:


#NEXUS

[Tree distribution from: The global diversity of birds in space and time;
W. Jetz, G. H. Thomas, J. B. Joy, K. Hartmann, A. O. Mooers
doi:10.1038/nature11631]
[Subsampled and pruned from birdtree.org on 2018-02-13T21:49:59Z ]
[Data: "Hackett All Species" (see Jetz et al. 2012 supplement for details)]

BEGIN TREES;
TREE tree_4104 =
TREE tree_9438 =
...
TREE tree_1868 =
TREE tree_0011 =
END;


Looks like your download did not work and the error message from read.nexus
could improve.

Cheers,
Klaus


On Tue, Feb 13, 2018 at 5:11 PM, Agus Camacho <agus.cama...@gmail.com>
wrote:

> Dear list,
>
> I am puzzled with this error when importing a bird tree subset from
> https://birdtree.org/subsets/ .
>
> The subset is supposed to be in nexus format. However,
>
> trees<-read.nexus("output.nex")
>
> Error in if (tp[3] != "") obj$node.label <- tp[3] :
>
> missing value where TRUE/FALSE needed
>
>
> read.tree gives the same error
> read.newick gets stuck, it never ends reading in the tree
>
> I tried scanning the file and then several notes appear before the actual
> trees Thus, I tried read.annotated.nexus from OutbreakTools.
> But it gives another error:
>
> Error in read.annotated.nexus("C:/Users/agusc/Dropbox/output.nex") :
> .annotated.clado.build is not yet implemented.
>
>
> I tried with several different subsets from different backbones and none
> worked.
> Here is one: https://www.dropbox.com/s/037u6pjr4006o6r/output.nex?dl=0
>
> Somebody has reported this error previously in the internet, but it still
> seems to be unsolved.
>
> Thanks in advance and best regards.
> Agus
>
>
> --
> Agustín Camacho Guerrero.
> http://www.agustincamacho.com
> Doutor em Zoologia.
> Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
> Biociências, USP.
> Rua do Matão, trav. 14, nº 321, Cidade Universitária,
> São Paulo - SP, CEP: 05508-090, Brasil.
>
> [[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-ph...@r-project.org/
>



-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston
http://www.phangorn.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/

Re: [R-sig-phylo] Rescaling a cophenetic matrix based on the Early Burst Model?

[Klaus Schliep]

Dear Max,
there is a function called coph in phangorn. The function is hidden in the
name space, so you have to use phangorn:::coph.
It returns a dist object instead of a matrix and computing time is roughly
O(n^2) instead of O(n^3).
> library(phangorn)> tree <- rtree(1)
> system.time(cm1 <- cophenetic(tree))
   user  system elapsed
  4.312   0.428   4.752
> system.time(cm2 <- phangorn:::coph(tree))
   user  system elapsed
  0.228   0.036   0.262
For large trees you are likely running out of memory, so Cecile's approach
is in the long run probably better.
Cheers,
Klaus



On Fri, Dec 15, 2017 at 5:26 PM, Cecile Ane <cecile@wisc.edu> wrote:

> Not in phylolm, because it’s faster to avoid manipulating this large
> matrix, if possible (and it is indeed possible for many purposes). If you
> can rephrase your calculations to use branching times (distance from a node
> to its descendant tips) or using the nodes’ distance to the root, then you
> could use functions like pruningwise.branching.times and
> pruningwise.distFromRoot in phylolm. Other folks could chime in, for tools
> in other packages.
> Cécile
>
> On Dec 15, 2017, at 3:55 PM, Max Farrell <maxwellfarr...@gmail.com maxwellfarr...@gmail.com>> wrote:
>
> Hi Cecile,
>
> Thanks for the input - it looks like this function does perform the tree
> rescaling, and in about 1/4 of the time, so this could help speed up the
> code!
>
> However, the function doesn't seem to return a covariance matrix. Is there
> a way to get the cophenetic matrix with the out put of this function?
> Otherwise I'm stuck using cophenetic again...
>
> Max
>
> Max
>
> On Fri, Dec 15, 2017 at 4:36 PM, Cecile Ane <cecile@wisc.edu<mailto:ce
> cile@wisc.edu>> wrote:
> In the package phylolm, the function "transf.branch.lengths” might do what
> you need. It has an option model=“EB” for early burst.
> Cécile
>
> > On Dec 15, 2017, at 3:27 PM, Max Farrell <maxwellfarr...@gmail.com<
> mailto:maxwellfarr...@gmail.com>> wrote:
> >
> > I have been using the rescale function from geiger for a link prediction
> > model I recently helped develop (https://arxiv.org/abs/1707.08354). I'm
> now
> > running this model on a much larger dataset and part of the code is
> > computing cophenetic(rescale(phy)) with 'EB' rescaling many many times.
> >
> > We're finding that calling cophenetic() is the rate limiting step, and if
> > we can avoid this function we expect to speed up our code by up to 4
> times.
> > This would be very useful as our simulation has been running for over 50
> > days now...
> >
> > I was wondering if there is a way to do EB transformations on a
> > phylogenetic distance matrix directly so as to avoid using the cophenetic
> > function?
> >
> > Any help or insights would be greatly appreciated!
> >
> >   [[alternative HTML version deleted]]
> >
> > ___
> > R-sig-phylo mailing list - R-sig-phylo@r-project.org 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-ph...@r-project.org/
>
> ___
> R-sig-phylo mailing list - R-sig-phylo@r-project.org 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-ph...@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-ph...@r-project.org/
>



-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston
http://www.phangorn.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/

Re: [R-sig-phylo] Possible Bug in ape::read.tree

[Klaus Schliep]

t;0:ID=AN45],(AN90:0.984,
>>> AN91:0.893):0.158[&:Ev=1>0:ID=AN89]):0.262[&:Ev=0>
>>> 1:S=Fungi:ID=AN44]):0.394[&:Ev=0>1:S=Opisthokonts:ID=
>>> AN3],(AN93:0.457,AN94:0.462):0.860[&:Ev=0>1:S=
>>> Dictyostelium:ID=AN92]):0.301[&:Ev=0>1:S=Unikonts:ID=
>>> AN2],AN95:1.440,AN96:2.000):2.000[&:Ev=0>1:S=Eukaryota:
>>> ID=AN1])[&:Ev=1>0:ID=AN0];
>>>
>>> I get the following error: "The tree has apparently singleton node(s):
>>> cannot read tree file.  Reading Newick file aborted at tree no. 1",
>>> which seems to be wrong since I can actually read this tree without
>>> problem
>>> using rncl::read_newick_phylo and this online tool
>>> http://etetoolkit.org/treeview/?treeid=6e192f20a3226cfdde219531cc533e
>>> 9f=ce443ed1e53858bf4e11d1e069c7a927
>>>
>>> I understand that this particular tree is a modified version of Newick's
>>> format, but this is the first time that I have problems reading this type
>>> of tree.
>>>
>>> Here my session info:
>>>
>>> devtools::session_info()
>>>>
>>> Session info
>>> 
>>> -
>>>  setting  value
>>>  version  R version 3.4.1 (2017-06-30)
>>>  system   x86_64, linux-gnu
>>>  ui   RStudio (1.0.143)
>>>  language (EN)
>>>  collate  en_US.UTF-8
>>>  tz   America/New_York
>>>  date 2017-09-05
>>>
>>> Packages
>>> 
>>> -
>>>  package * version date   source
>>>  ape   4.1 2017-02-14 CRAN (R 3.4.0)
>>>  assertthat0.2.0   2017-04-11 CRAN (R 3.4.0)
>>>  base* 3.4.1   2017-06-30 local
>>>  compiler  3.4.1   2017-06-30 local
>>>  datasets* 3.4.1   2017-06-30 local
>>>  devtools  1.13.3  2017-08-02 CRAN (R 3.4.0)
>>>  digest0.6.12  2017-01-27 CRAN (R 3.4.0)
>>>  graphics* 3.4.1   2017-06-30 local
>>>  grDevices   * 3.4.1   2017-06-30 local
>>>  grid  3.4.1   2017-06-30 local
>>>  lattice   0.20-35 2017-03-25 CRAN (R 3.4.1)
>>>  magrittr  1.5 2014-11-22 CRAN (R 3.4.0)
>>>  memoise   1.1.0   2017-04-21 CRAN (R 3.4.0)
>>>  methods * 3.4.1   2017-06-30 local
>>>  nlme  3.1-131 2017-02-06 CRAN (R 3.4.1)
>>>  parallel  3.4.1   2017-06-30 local
>>>  prettyunits   1.0.2   2015-07-13 CRAN (R 3.4.1)
>>>  progress  1.1.2   2016-12-14 CRAN (R 3.4.1)
>>>  R62.2.2   2017-06-17 CRAN (R 3.4.0)
>>>  Rcpp  0.12.12 2017-07-15 CRAN (R 3.4.0)
>>>  rncl  0.8.2   2016-12-16 CRAN (R 3.4.1)
>>>  stats   * 3.4.1   2017-06-30 local
>>>  tools 3.4.1   2017-06-30 local
>>>  utils   * 3.4.1   2017-06-30 local
>>>  withr 2.0.0   2017-07-28 CRAN (R 3.4.0)
>>>
>>>
>>> Any ideas?
>>>
>>> Best,
>>>
>>> George G. Vega Yon
>>> +1 (626) 381 8171
>>> http://cana.usc.edu/vegayon
>>>
>>> [[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-ph...@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-ph...@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-ph...@r-project.org/
>>
>>
> ___
> 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-ph...@r-project.org/
>



-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston
http://www.phangorn.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/

Re: [R-sig-phylo] selecting a set of incongruent trees from a posterior distribution

[Klaus Schliep]

gt; > > haven't found anything.
> > >
> > > I running mixed models, which take a lot of computational space on my
> lap
> > > top. In effort to account for phylogenetic uncertainty, without having
> to
> > > run 100s of chains, I figured maybe i could run models across a (much)
> > > shorter list that accounts for more diversity in tree shape observed
> > within
> > > the posterior distribution. Not sure if this makes sense, and/or is
> > > extremely complicated?
> > >
> > > Thanks for you time!
> > >
> > > Jesse
> > >
> > > [[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-ph...@r-project.org/
> > >
> >
> >     [[alternative HTML version deleted]]
> >
> >
> >
> > --
> >
> > Subject: Digest Footer
> >
> > ___
> > R-sig-phylo mailing list
> > R-sig-phylo@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
> >
> > --
> >
> > End of R-sig-phylo Digest, Vol 114, Issue 10
> > 
> >
> --
> ==
> Santiago Sanchez-Ramirez, PhD
> Postdoctoral Associate
> Ecology and Evolutionary Biology
> University of Toronto
> ==
>
> [[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-ph...@r-project.org/
>



-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston
http://www.phangorn.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/

Re: [R-sig-phylo] A possible alternate MRCA function to APE's getMRCA

[Klaus Schliep]

gt;
>>
>> ## traverse lineages for remaining tips, stop if hit common ancestor
>> for (i in 2:length(tnd)) {
>> done <- FALSE
>> cnd <- tnd[i]
>> while (!done) {
>> cpar<- pvec[cnd] # get immediate parent
>> if (cpar %in% pars) {
>> if (cpar == rootnd) return(rootnd) # early exit
>> mrcand <- min(mrcand,cpar); # shallowest node has lowest nodeid
>>
>> done <- TRUE
>> }
>> cnd <- cpar # keep going!
>> }
>> }
>> mrcand
>> }
>>
>> I also played around with compiling the function:
>>
>> library(compiler);
>> getMRCA_em_cmp<- cmpfun(getMRCA_em);
>> getMRCA_new_cmp<- cmpfun(getMRCA_new);
>>
>> Here are some timings (*_em is your code, *_new is the code above). For
>> my original example (100K tips, 500 tips in MRCA query) there is not a
>> whole lot of difference:
>>
>> microbenchmark(getMRCA_em(phy, tips), getMRCA_new(phy, tips),
>> getMRCA_em_cmp(phy, tips), getMRCA_new_cmp(phy, tips));
>> Unit: milliseconds
>> expr min lq mean median uq max neval
>>getMRCA_em(phy, tips) 23.9987929.8196433.6743732.391
>> 4935.61888136.34521100
>>   getMRCA_new(phy, tips) 23.7694128.0250031.8428531.370
>> 0134.8964044.66380100
>>getMRCA_em_cmp(phy, tips) 20.8192423.6011828.2395728.252
>> 0230.1853559.72102100
>>   getMRCA_new_cmp(phy, tips) 20.1658823.1957127.1961026.657
>> 8329.8214942.61702100
>>
>> However, for the worst-case scenario that Klaus devised the results are
>> very impressive:
>>
>> microbenchmark(getMRCA_em(tree, c(1,2)), getMRCA_new(tree, c(1,2)),
>> getMRCA_em_cmp(tree, c(1,2)), getMRCA_new_cmp(tree, c(1,2)));
>> Unit: milliseconds
>> expr min lq mean median uq max neval
>>getMRCA_em(tree, c(1, 2)) 124.79040140.44473146.98292145
>> .39224150.96044286.78335100
>>   getMRCA_new(tree, c(1, 2)) 123.39924139.41609145.21311143
>> .92062148.34205297.73846100
>>getMRCA_em_cmp(tree, c(1, 2)) 23.7210526.8494429.9495029.879
>> 3532.2649041.99264100
>>   getMRCA_new_cmp(tree, c(1, 2)) 23.1522325.7468629.2152629.253
>> 7030.8091569.68522100
>>
>> So it might be worth compiling.
>>
>> 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 10 Jun, 2017, at 05:26, Emmanuel Paradis <emmanuel.para...@ird.fr
>>> <mailto:emmanuel.para...@ird.fr>> wrote:
>>>
>>> Joseph, Klaus,
>>>
>>> This is great. I slightly edited the code and renamed the argument
>>> 'tips' to 'tip' (as in the original definition). I run on a bunch of random
>>> trees and the results are exactly identical with the current version of
>>> getMRCA. I'm changing the code in ape now.
>>>
>>> Cheers,
>>>
>>> Emmanuel
>>>
>>> getMRCA <- function(phy, tip)
>>> {
>>>if (!inherits(phy, "phylo"))
>>>stop('object "phy" is not of class "phylo"')
>>>if (length(tip) < 2) return(NULL)
>>>Ntip <- length(phy$tip.label)
>>>rootnd <- Ntip + 1L
>>>
>>>pars <- integer(phy$Nnode) # worst case assignment, usually far too
>>> long
>>>mrcaind <- 1L # index in pars of the mrca node. will return highest
>>> one traversed by other lineages
>>>tnd <- if (is.character(tip)) match(tip, phy$tip.label) else tip
>>>
>>>## build a lookup table to get parents faster
>>>pvec <- integer(Ntip + phy$Nnode)
>>>pvec[phy$edge[, 2]] <- phy$edge[, 1]
>>>
>>>## get entire lineage for first tip
>>>nd <- tnd[1]
>>>for (k in 1:phy$Nnode) {
>>>nd <- pvec[nd]
>>>pars[k] <- nd
>>>if (nd == rootnd) break
>>>}
>>>pars <- pars[1:k] # delete the rest
>>>
>>>## traverse lineages for remaining tips, stop if hit common ancestor
>>>for (i in 2:length(tnd)) {
>>>done <- FALSE
>>>cnd <- tnd[i]
>>>while (!done) {
>>>cpar <- pvec[cnd] # get immediate parent
>>>if (cpar %in%

Re: [R-sig-phylo] A possible alternate MRCA function to APE\'s getMRCA

[Klaus Schliep]

our
> >>> function.
> >>>
> >>> Cheers
> >>>
> >>> Juan
> >>>
> >>> --
> >>>
> >>> Dr. Juan A. Balbuena
> >>> Cavanilles Institute of Biodiversity and Evolutionary Biology
> >>> University of Valencia
> >>> http://www.uv.es/~balbuena <http://www.uv.es/%7Ebalbuena>
> >>> P.O. Box 22085
> >>> http://www.uv.es/cophylpaco
> >>> <http://www.uv.es/cavanilles/zoomarin/index.htm>
> >>> 46071 Valencia, Spain
> >>> e-mail: j.a.balbu...@uv.es <mailto:j.a.balbu...@uv.es>tel. +34 963
> >>> 543 658fax +34 963 543 733
> >>> 
> >>> *NOTE!*For shipments by EXPRESS COURIER use the following street
> address:
> >>> C/ Catedrático José Beltrán 2, 46980 Paterna (Valencia), Spain.
> >>> 
> >>>
> >>>
> >>> <https://www.avast.com/sig-email?utm_medium=email_
> source=link_campaign=sig-email_content=emailclient>
> >>>
> >>>  Libre de virus. www.avast.com
> >>> <https://www.avast.com/sig-email?utm_medium=email_
> source=link_campaign=sig-email_content=emailclient>
> >>>
> >>>
> >>>
> >>> <#DAB4FAD8-2DD7-40BB-A1B8-4E2AA1F9FDF2>
> >>>
> >>>
> >>> ___
> >>> 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/
> >>>
> >>
> >> ___
> >> 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-ph...@r-project.org/
>



-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston
http://www.phangorn.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/

Re: [R-sig-phylo] A possible alternate MRCA function to APE's getMRCA

[Klaus Schliep]

p example
> set.seed(13);
> phy <- rtree(10);
> node <- 13; # make sure MRCA is not the root
> tt <- unlist(Descendants(phy, 13));
> set.seed(13);
> tips <- sample(tt, 50);
>
> # APE
> system.time(nd.ape <- getMRCA(phy, tips)); nd.ape;
>user  system elapsed
>  14.459   0.066  14.482
> [1] 13
>
> # ALT
> system.time(nd.alt <- alt_mrca(phy, tips)); nd.alt;
>user  system elapsed
>   2.557   0.534   3.063
> [1] 13
>
> Ok, that was for a shallow MRCA node (i.e. close to the root). How does it
> fare for recent MRCAs? Here is a timing for sister tips:
>
> # APE
> system.time(nd.ape <- getMRCA(phy, c(13,14))); nd.ape;
>user  system elapsed
>  14.715   0.063  14.727
> [1] 100027
>
> # ALT
> system.time(nd.alt <- alt_mrca(phy, c(13,14))); nd.alt;
>user  system elapsed
>   0.047   0.015   0.061
> [1] 100027
>
> Way faster! Getting the entire lineage of the first tip does not seem to
> hurt us.
>
> We can also look how the timing scales as the size of the query set
> changes (I've made sure that the MRCA is identical in every case):
> (APE is blue, alt_mrca is red). So the APE version seems insensitive to
> the number of tips query. I suppose there must be a point where the timings
> cross and the alt_mrca function takes more time, but I have not found it.
> Besides, if the query set becomes large enough then the MRCA is likely the
> root itself, and in that case we may exit early:
>
> # APE with just 2 tips
> system.time(nd.ape <- getMRCA(phy, c(1,10))); nd.ape;
>user  system elapsed
>  14.627   0.055  14.642
> [1] 11
>
> # ALT with a random 500 tips
> system.time(nd.alt <- alt_mrca(phy, sample(1:10, 500))); nd.alt;
>user  system elapsed
>   0.292   0.058   0.346
> [1] 11
>
> Not bad. The code for the alternate MRCA is below. I've considered
> lapply-ing this (i.e. across query nodes), but that would negate the
> potential early exit (as *all* tip-to-root lineages would have to be
> traversed). Note that in all of these timings above the APE getMRCA
> function uses C code but the alternative is pure R. I imagine the alt_mrca
> could be made even speedier if coded in C.
>
> HTH.
> Joseph.
>
> alt_mrca <- function (phy, tips) {
> rootnd <- length(phy$tip.label) + 1;
> pars <- numeric(); # parents in lineage of first queried tip
> mrcaind <- 1;  # index in pars of the mrca node
> tnd <- NULL;   # node ids of tips
> if ("character" %in% class(tips)) {
> tnd <- match(tips, phy$tip.label);
> } else {
> tnd <- tips;
> }
>
>
> # get entire lineage for first tip
> nd <- tnd[1];
> repeat {
> nd <- phy$edge[which(phy$edge[,2] == nd),1]; # get immediate
> parent
> pars <- c(pars, nd);
> if (nd == rootnd) {
> break;
> }
> }
>
>
> # traverse lineages for remaining tips, stop if hit common ancestor
> for (i in 2:length(tnd)) {
> done <- FALSE;
> cnd <- tnd[i];
> while (!done) {
> cpar <- phy$edge[which(phy$edge[,2] == cnd),1]; # get
> immediate parent
> if (cpar %in% pars) {
> if (cpar == rootnd) {
> # early exit! if the mrca of any set of taxa is the
> root then we are done
> return(rootnd);
> }
> idx <- which(pars == cpar);
> if (idx > mrcaind) {
> mrcaind <- idx;
> }
> done <- TRUE;
> } else {
> # keep going!
> cnd <- cpar;
> }
> }
> }
> return(pars[mrcaind]);
> }
>
> 
> 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
>
>
>
>
> ___
> 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-ph...@r-project.org/
>



-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston
http://www.phangorn.org/
___
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/

Re: [R-sig-phylo] Graphically comparing multiple trees

[Klaus Schliep]

Hi Vojtěch,


On Wed, Apr 26, 2017 at 1:09 PM, Vojtěch Zeisek <vo...@trapa.cz> wrote:

> Hi,
> thank You for the advice. I tried:
> densiTree(x=hirta, type="phylogram", scaleX=TRUE, width=2,
> col=rainbow(length(hirta)), font=3, underscore=FALSE)
> looks relatively nicely. I just wonder why last 2 options have no effect
> (i.e.
> I wish tip labels in italics and without underscore).


you are right the last 2 options have no effect, at least in the CRAN
version,
I added them and a few more on the in the development version on github:
devtools::install_github("KlausVigo/phangorn")



> I also didn't find how
> to add little vertical offset among trees to show identical lines in
> parallel.
>

I thought this would be your little homework, but Liam spoiled it ;) and
did all the work for you.
The densiTree function is just about 40 lines, so it would be easy just add
some offset, etc.
and I always look forward to pull requests.

Cheers,
Klaus


> Without it, the display is very messy and unreadable... :-(
> Sincerely,
> V.
>
> Dne úterý 25. dubna 2017 23:04:01 CEST jste napsal(a):
> > Hi Vojtěch,
> > you could try adopt densiTree() in phangorn to do this. It just plots
> each
> > tree separate over each other. You probably want to of move each tree a
> bit
> > up or down and add some more control for colors. It is based on the
> > internal function of from plot.phylo().
> > Cheers,
> > Klaus
> >
> > On Tue, Apr 25, 2017 at 2:25 PM, Liam J. Revell <liam.rev...@umb.edu>
> wrote:
> > > Hi Vojtěch.
> > >
> > > This kind of plot is possible to make using plotTree in phytools via
> the
> > > arguments tips (which sets the vertical coordinates of the tips,
> > > regardless
> > > of topology) & add (a logical value indicating whether to add the tree
> > > tree
> > > to an existing plot or to open a new plot). (All arguments for plotTree
> > > are
> > > documented in the man page for plotSimmap which plotTree uses
> internally.)
> > > Unfortunately, there is some complication. For instance, unless we set
> > > xlim
> > > & ylim, then plotTree will reset these values for each new plotted
> tree.
> > >
> > > Later this evening I will do my best to make an example & post it to my
> > > blog.
> > >
> > > - Liam
> > >
> > > On 4/25/2017 12:20 PM, Vojtěch Zeisek wrote:
> > >> Hello,
> > >> for comparison of two trees I can use very nice function cophyloplot
> > >> plotting
> > >> two trees (left and right) and connecting respective nodes by lines.
> Very
> > >> nice
> > >> and convenient to read. But only for two trees. Displaying multiple
> trees
> > >> in multiple comparisons is not very convenient.
> > >> To display dozens to hundreds of trees there is densitree. Also nice,
> but
> > >> for this purpose I don't like its display.
> > >> I have several trees (~5) and I wish to compare their topologies, show
> > >> supports (at least for differing nodes) and highlight differences. I
> > >> thought
> > >> about some overlay/parallel plotting (similar to the attached image)
> > >> where
> > >> there would be complete topologies displayed and incongruences would
> be
> > >> easily
> > >> visible. It would be probably doable by plotting all separate trees by
> > >> plot.phylo and then combining and tuning the figure in some vector
> editor
> > >> (like Inkscape). But I hope there is some more automated way to do it.
> > >> :-)
> > >> Sincerely,
> > >> V.
> --
> Vojtěch Zeisek
> https://trapa.cz/en/
>
> Department of Botany, Faculty of Science
> Charles University, Prague, Czech Republic
> https://www.natur.cuni.cz/biology/botany/
>
> Institute of Botany, Czech Academy of Sciences
> Průhonice, Czech Republic
> http://www.ibot.cas.cz/en/
>
> ___
> 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-ph...@r-project.org/
>



-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston
http://www.phangorn.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/

Re: [R-sig-phylo] Graphically comparing multiple trees

[Klaus Schliep]

Hi Vojtěch,
you could try adopt densiTree() in phangorn to do this. It just plots each
tree separate over each other. You probably want to of move each tree a bit
up or down and add some more control for colors. It is based on the
internal function of from plot.phylo().
Cheers,
Klaus

On Tue, Apr 25, 2017 at 2:25 PM, Liam J. Revell <liam.rev...@umb.edu> wrote:

> Hi Vojtěch.
>
> This kind of plot is possible to make using plotTree in phytools via the
> arguments tips (which sets the vertical coordinates of the tips, regardless
> of topology) & add (a logical value indicating whether to add the tree tree
> to an existing plot or to open a new plot). (All arguments for plotTree are
> documented in the man page for plotSimmap which plotTree uses internally.)
> Unfortunately, there is some complication. For instance, unless we set xlim
> & ylim, then plotTree will reset these values for each new plotted tree.
>
> Later this evening I will do my best to make an example & post it to my
> blog.
>
> - Liam
>
> Liam J. Revell, Associate Professor of Biology
> University of Massachusetts Boston
> web: http://faculty.umb.edu/liam.revell/
> email: liam.rev...@umb.edu
> blog: http://blog.phytools.org
>
> On 4/25/2017 12:20 PM, Vojtěch Zeisek wrote:
>
>> Hello,
>> for comparison of two trees I can use very nice function cophyloplot
>> plotting
>> two trees (left and right) and connecting respective nodes by lines. Very
>> nice
>> and convenient to read. But only for two trees. Displaying multiple trees
>> in
>> multiple comparisons is not very convenient.
>> To display dozens to hundreds of trees there is densitree. Also nice, but
>> for
>> this purpose I don't like its display.
>> I have several trees (~5) and I wish to compare their topologies, show
>> supports (at least for differing nodes) and highlight differences. I
>> thought
>> about some overlay/parallel plotting (similar to the attached image) where
>> there would be complete topologies displayed and incongruences would be
>> easily
>> visible. It would be probably doable by plotting all separate trees by
>> plot.phylo and then combining and tuning the figure in some vector editor
>> (like Inkscape). But I hope there is some more automated way to do it. :-)
>> Sincerely,
>> V.
>>
>>
>>
>> ___
>> 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-ph...@r-project.org/
>>
>>
> ___
> 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-ph...@r-project.org/
>



-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston
http://www.phangorn.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/

Re: [R-sig-phylo] distances between the nodes

[Klaus Schliep]

Dear Justyna

tr$edge does not store coordinates, see help(phylo).
if you want the coordinates of the nodes of your last plotted tree you can
get them:
library(ape)
tree = rtree(5)
nodelabels()
tiplabels()
last_plot = get("last_plot.phylo", envir = .PlotPhyloEnv)
last_plot$xx and last_plot$yy contain the coordinates of the nodes.

there is also a (generic) function cophenetic() in ape to compute the
patristic distance and for larger trees this is much faster than distTips.

Regards,
Klaus


On Fri, Feb 24, 2017 at 11:03 AM, Eduardo Ascarrunz <ear...@gmail.com>
wrote:

> Hi Justyna,
>
> If you want patristic distances (or node distances) you can use distTips
> from the adephylo package.
>
> To make edgelabels show the branch lengths you can use
> edgelabels(tr$edge.length)
>
> Does that do what you wanted?
>
> Cheers,
>
> Eduardo
>
> 2017-02-23 16:41 GMT+01:00 Wierzbinska, Justyna <j.wierzbinska@dkfz-
> heidelberg.de>:
>
>>
>>
>> Dear Ape users,
>>
>>
>>
>> I’be grateful for your help with my ape problems.
>>
>> I find it difficult to figure out how I can look at the distances between
>> the main nodes generated in the phylogenetic tree below. My idea is to
>> extract the distances between NBCs to hiMBCs. Somehow looking at the
>> exemplary object tr and using tr$edge and tr$edge.length it doesn’t seem to
>> be an easy task.
>>
>>
>>
>>
>>
>> *[image: cid:part1.C97E7A27.01EB5ADD@uq.edu.au]*
>>
>>
>>
>> I’ve tried to add node labels using nodelabels() and edge labels but
>> somehow the assigned labels are different then the coordinates stored in
>> the tr$edge would suggest (when I look at the tip projection).
>>
>> Any ideas how I can do it in an automatic way? Is it just a plotting
>> problem that the coordinates of edges don’t correspond to what I see in the
>> tree.
>>
>>
>>
>>
>>
>> Thank you very much.
>>
>> Kind regards,
>>
>> Justyna
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> Division of Epigenomics and Cancer Risk Factors (C010)
>>
>>
>>
>> *German Cancer Research Center (DKFZ)*
>>
>> Foundation under Public Law
>>
>> Im Neuenheimer Feld 280
>>
>> 69120 Heidelberg
>>
>> Germany
>>
>> +49 6221 42 <+49%206221%2042> 3359
>>
>> j.wierzbin...@dkfz-heidelberg.de
>>
>>
>>
>> Management Board: Prof. Dr. Michael Baumann, Prof. Dr. Josef Puchta
>>
>> VAT-ID No.: DE143293537
>>
>>
>>
>> ___
>> 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-ph...@r-project.org/
>>
>
>
> ___
> 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-ph...@r-project.org/
>



-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston
http://www.phangorn.org/
___
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/

Re: [R-sig-phylo] Enigmatic error with drop.tip

[Klaus Schliep]

Hi Juan,
it seems that you trim too many taxa.

I can replicate the error message if I trim for example 19 out of a 20 tips
of a tree.
drop.tip(tree, sample(20, 19))
Error in integer(max(oldnodes)) : vector size cannot be infinite
In addition: Warning message:
In max(oldnodes) : no non-missing arguments to max; returning -Inf
> traceback()
3: integer(max(oldnodes))
2: collapse.singles(phy)
1: drop.tip(tree, sample(20, 19))

Maybe you include a check into your function z that the pruned trees have
at least two tips left
if( length( setdiff(treeP$tip.label, colnames(hp))) < ( length(
treeP$tip.label ) -2 ) )

Cheers,
Klaus



On Fri, Jan 20, 2017 at 7:26 AM, Juan Antonio Balbuena <j.a.balbu...@uv.es>
wrote:

> Hello all
>
> I am working with a cophylogenetic problem where I have a host tree, a
> parasite tree and a binary matrix HP encapsulating the host-parasite
> associations . After trimming the matrix according to a given criterion, I
> require to prune the trees accordingly. Since I am getting errors at the
> execution I devised the following control function to check the number of
> tree nodes after pruning:
> z <- function (hp, treeH,treeP) {
>   treeh <- drop.tip(treeH, setdiff(treeH$tip.label, rownames(hp)))
>   treep <- drop.tip(treeP, setdiff(treeP$tip.label, colnames(hp)))
>   res <- c(treeh$Nnode,treep$Nnode)
>  return(res)
>  }
>
> Then I try the following
>
> x <- sapply (THP, z, treeH=TreeH[[7]], treeP= TreeP[[7]])
>
> where THP is a list of trimmed HP matrices and TreeH and TreeP are
> multiphylo objects and get the following error:
>
> Error in integer(max(oldnodes)) : vector size cannot be infinite
>
> The traceback is
>
> 6. integer(max(oldnodes))
> 5. collapse.singles(phy)
> 4. drop.tip(treeP, setdiff(treeP$tip.label, colnames(hp)))
> 3. FUN(X[[i]], ...)
> 2. lapply(X = X, FUN = FUN, ...)
> 1. sapply(THP, z, treeH = TreeH[[7]], treeP = TreeP[[7]])
>
> This error appears only with this particular pair of host and parasite
> trees (#7) and I cannot find a reason why this case should be different.
>
> Any help will be very much appreciated.
>
> Juan A. Balbuena
>
> --
>
> Dr. Juan A. Balbuena
> Cavanilles Institute of Biodiversity and Evolutionary Biology
> University of Valencia
> http://www.uv.es/~balbuena
> P.O. Box 22085
> http://www.uv.es/cophylpaco
> <http://www.uv.es/cavanilles/zoomarin/index.htm>
> 46071 Valencia, Spain
> e-mail: j.a.balbu...@uv.estel. +34 963 543 658
> <+34%20963%2054%2036%2058>fax +34 963 543 733
> <+34%20963%2054%2037%2033>
> 
> *NOTE!* For shipments by EXPRESS COURIER use the following street address:
> C/ Catedrático José Beltrán 2, 46980 Paterna (Valencia), Spain.
> 
>
> _______
> 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-ph...@r-project.org/
>



-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston
http://www.phangorn.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/

Re: [R-sig-phylo] multi2di fails on multiPhylo objects

[Klaus Schliep]

Hi Yan,

it seems we introduced a small bug with making these functions generic.
Replace in the function .multi2di_ape the line
phy <- reorder(phy)
with
phy <- .reorder_ape(phy, "cladewise", FALSE, n, 1L)
and it should work. You have to recompile ape afterwards.

Cheers,
Klaus





On Fri, Jan 13, 2017 at 1:24 PM, Yan Wong <y...@yanwong.me> wrote:

> Thanks. I think that’s because read_nexus_phylo() doesn’t do
> .compressTipLabel automatically
>
> n <- read_nexus_phylo("tmp.nex")
> n2 <- .compressTipLabel(n)
> multi2di(n2) #this causes an error.
>
> Unfortunately I need the .compressTipLabel functionality really.
>
> Yan
>
>
> On 13 Jan 2017, at 17:06, François Michonneau <
> francois.michonn...@gmail.com> wrote:
>
> > Hi Yan,
> >
> >  I'm not sure about what's going on here, but it works with the rncl
> > package (http://cran.r-project.org/web/packages/rncl/index.html)
> >
> >  install.packages("rncl")
> >  library(rncl)
> >  n <- read_nexus_phylo("tmp.nex")
>
> ___
> 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-ph...@r-project.org/
>



-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston
http://www.phangorn.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/

Re: [R-sig-phylo] phylo from igraph

[Klaus Schliep]

Hi Giulio,
there is an as_phylo function in igraph and as.igraph in ape.
Cheers,
Klaus




On Fri, Jan 6, 2017 at 5:02 PM, Giulio V. Dalla Riva <gv...@uclive.ac.nz>
wrote:

> Dear Phyloers,
>
>
> Has anybody already implemented a function to convert an igraph graph
> object into a phylo object?
>
>
> (Let's pretend we are in the best of the worlds and the graph object is
> indeed a tree and all).
>
>
> Best,
>
>
> Giulio Valentino Dalla Riva
>
> Beaty Biodiversity Research Centre
> University of British Columbia
> Vancouver, Canada
>
> [[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-ph...@r-project.org/
>



-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston
http://www.phangorn.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/

Re: [R-sig-phylo] Midpoint root a tree with root() function ape ver 4.0

[Klaus Schliep]

Hi Knud,
in phangorn I try to take care that the node labels are assigned to the
right nodes after using midpoint, see attached code and pic. It seems to
work quite nicely. It would be useful if you could supply a reproducible
example or the tree you have problems with and open otherwise an issue on
github.
Regards,
Klaus




On Thu, Dec 15, 2016 at 9:14 AM, Todd Knutson <knut0...@umn.edu> wrote:

> Hi,
>
> Is there any way to use the updated root() function in ape ver 4.0 to find
> and set the midpoint root in a tree? I greatly appreciate the updated
> root() function to include the “edgelabel = TRUE” option, so that when I
> have bootstrapping support values listed as node labels, they get assigned
> to the proper edge of the tree after re-rooting.
>
> However, using the midpoint.root() function from phytools or midpoint()
> from the phangorn packages alter the node labels after rooting, and the
> bootstrapping values get assigned to the wrong edges.
>
> Thus, I would love to use the new root() function, with “edgelabel = TRUE”
> option enabled, to midpoint a tree. I would appreciate any suggestions.
>
> Thanks,
> Todd
> ___
> 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-ph...@r-project.org/




-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston
http://www.phangorn.org/
library(phangorn)
library(ape)

data("Laurasiatherian")
tree <- nj(dist.ml(Laurasiatherian))
set.seed(42)
trees <- bootstrap.phyDat(Laurasiatherian, function(x)nj(dist.ml(x)))

tree1 <- plotBS(tree, trees, "phylogram")
pdf("Bootstrap.pdf")
par(mar=c(1,1,1,1))
par(mfrow=c(3,1))
plot(tree1, show.node.label = TRUE)
tree2 <- midpoint(tree1)
plot(tree2, show.node.label = TRUE)
tree$tip.label[43]
tree3 = tree1
# make edge to Cat large
ind = which(tree1$edge[,2]==43)
tree3$edge.length[ind]=.5
tree3 = midpoint(tree3)
plot(tree3, show.node.label = TRUE)
dev.off()


Bootstrap.pdf
Description: Adobe PDF document
___
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/

Re: [R-sig-phylo] compressTipLabel as an option to read.trees()

[Klaus Schliep]

Hi Yan,
Joseph was right. In read.nexus you need a TRANSLATE block, just a
TAXLABELS is not enough. Then read.nexus returns the compressed object and
is 10x faster to read in (for 1000 trees with 1000 taxa on my machine).
There is also the package rncl (Nexus Class Library), it is faster to read
in, even the pure R implementation with the TRANSLATE block is almost as
fast.
However the objects are actually quite a bit larger. It also stores the
edge matrix as doubles, and which I find dangerous.
Cheers,
Klaus

On Wed, Dec 14, 2016 at 4:44 PM, Yan Wong <y...@yanwong.me> wrote:

>
> On 14 Dec 2016, at 20:57, Emmanuel Paradis <emmanuel.para...@ird.fr>
> wrote:
>
> > What is the size of your problem?
>
> Erm, quite large. I am looking at tree comparison metrics for roughly
> 10,000 trees with perhaps 10,000 tips on each, replicated several times.
> The newick files themselves take up gigabyes uncompressed. For this sized
> problem I’m likely to implement my own comparison metrics, but I want to
> trial this out with a tested library before rolling my own.
>
> > Do you use a recent version of ape? This function was improved one or
> two years ago.
>
> Yes, 4.0.
>
> But I’m happy for the moment to just leave this stuff running for days on
> a server, so it was just a quick suggestion really.
>
> Thanks for the quick reply
>
> Yan
> ___
> 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-ph...@r-project.org/
>



-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston
http://www.phangorn.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/

Re: [R-sig-phylo] concatenating loci with different taxa

[Klaus Schliep]

Dear Karla,

what is exactly the problem and what did you try? There is a (generic)
function cbind for "phyDat" objects in phangorn, which should exactly do
this task.

library(phangorn)
data(Laurasiatherian)
(x1 <- subset(Laurasiatherian, 1:20))
(x2 <- subset(Laurasiatherian, 11:30))
cbind(x1, x2)

As Emmanuel mentioned the package apex simplifies the whole process.

library(apex)
files <- dir(system.file(package="apex"),patter="patr", full=TRUE)
files
dat <- read.multiphyDat(files, format="FASTA")
concatenate(dat)

Regards,
Klaus




On Fri, Nov 18, 2016 at 6:21 AM, Karla Shikev <karlashi...@gmail.com> wrote:

> Dear all,
>
> I got alignments for many loci as nexus files and I've been trying to
> concatenate them into a single alignment. However, some taxa are missing
> from some of the loci, so that simple alternatives such as cbinding
> individual files in phangorn doesn't work.
>
> Any help with this will be greatly appreciated.
>
> Karla
>
> [[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-ph...@r-project.org/
>



-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston
http://www.phangorn.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/

Re: [R-sig-phylo] Anyone knows how to concatenate aligned genes sequences so as to create whole genome alignments?

[Klaus Schliep]

Dear Rav,
write.dna() from ape just does this.
Klaus

On Sep 12, 2016 9:07 AM, "Bhuller, Ravneet" <
ravneet.bhulle...@imperial.ac.uk> wrote:

> Hi Thibaut,
>
> Is there anyway I can save the concatenated alignment (created using apex
> and the concatenate function) in FASTA format?
>
> Regards,
>
> Rav
>
>
> On 12 Sep 2016, at 13:45, Thibaut Jombart  mailto:thibautjomb...@gmail.com>> wrote:
>
> Hi there,
>
> apex can do this using the 'concatenate' function:
> https://github.com/thibautjombart/apex
>
> Cheers
> Thibaut
>
>
> --
> Dr Thibaut Jombart
> Lecturer, Department of Infectious Disease Epidemiology, Imperial College
> London
> Head of RECON: https://reconhub.github.io/
> https://sites.google.com/site/thibautjombart/
> https://github.com/thibautjombart
> Twitter: @TeebzR
>
> On 12 September 2016 at 11:41, Bhuller, Ravneet <
> ravneet.bhulle...@imperial.ac.uk>
> wrote:
> Dear Members,
>
> Any suggestions on how to concatenate the aligned gene sequences in fasta
> format so as to get whole genome alignments?
> I need whole genome alignments as an input to a phylogenetic tool.
>
> Many thanks,
>
> Rav
>
> ___
> R-sig-phylo mailing list - R-sig-phylo@r-project.org 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-ph...@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-ph...@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/

Re: [R-sig-phylo] Ultrametric tree not recognized

[Klaus Schliep]

Hi Emmanuel,

it wasn't so much about the precision, but that it is not invariant against
linear transformation of the edge length.
At the moment the criterion inside is.ultrametric is the following:
var(xx[1:n])
where xx[1:n] is a vector with the tip-to-toot differences. And
is.ultrametric returns FALSE if the result of the expression larger than
tol.
What I suggest is doing something like this:
( max(xx[1:n]) - min(xx[1:n]) )  / max(xx[1:n])

Here is a small example:
set.seed(42)
tree <- read.tree(text=write.tree(rcoal(100), digits=4))
is.ultrametric(tree) # returns FALSE

And the 2 criteria are:
var(xx[1:n])
[1] 2.078812e-08
( max(xx[1:n]) - min(xx[1:n]) )  / max(xx[1:n])
[1] 0.0002541798

tree2 <- tree
tree2$edge.length[] = tree2$edge.length * 1000
is.ultrametric(tree2)

var(xx[1:n])
[1] 0.02078812
( max(xx[1:n]) - min(xx[1:n]) )  / max(xx[1:n])
[1] 0.0002541798

This is than at least invariant against transformation of the edge length.
You can possibly come up with other invariant criteria which may are even
more appropriate for this case.

Cheers,
Klaus





On Wed, Aug 17, 2016 at 3:04 AM, Emmanuel Paradis <emmanuel.para...@ird.fr>
wrote:

> Hi Klaus,
>
> The default is: is.ultrametric(phy, tol = .Machine$double.eps^0.5).
> According to ?.Machine this value does not depend on the machine precision.
> On my laptop, it gives:
>
> R> .Machine$double.eps^0.5
> [1] 1.490116e-08
>
> Best,
>
> Emmanuel
>
>
> Le 16/08/2016 à 18:15, Klaus Schliep a écrit :
>
>> We should consider to improve is.ultrametric to use as tolerance a certain
>> number significant figures (e.g. 8-10 digits). This way tol would only
>> depend on the tip-to-root distance and not a machine dependent minimal
>> value. I find the use of significant figures suits really well as a
>> stopping criterion for example in my ML functions, especially as it works
>> for small and large trees equally well.
>>
>> write.tree has an option to define the number of digits, it will be
>> generally less than the internal representation and I suspect one could
>> possible export with higher precision from other programs, but will result
>> in larger file sizes which one has to put in consideration for MCMC /
>> bootstrap tree files.
>> As a side note on the developer site of the r-project.org homepage (not
>> CRAN) is a link to this nice ancient paper (
>> http://www.validlab.com/goldberg/paper.pdf)
>>
>> Regards,
>> Klaus
>>
>>
>> PS: And Martin of course I am happy with Liam's solution as it uses
>> phangorn ;)
>>
>> On Tue, Aug 16, 2016 at 10:22 AM, Joseph W. Brown <josep...@umich.edu>
>> wrote:
>>
>> That’s a good tip, but the issue here is that many packages use
>>> is.ultrametric internally with a hard-coded (likely, default) value of
>>> tol.
>>> In this case, Martin simply cannot perform the analyses he wants to run
>>> because this option is inaccessible.
>>>
>>> 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 16 Aug, 2016, at 10:15, Klaus Schliep <klaus.schl...@gmail.com>
>>> wrote:
>>>
>>> Hello all,
>>> this may come be surprising to many, but consulting the manual
>>> ?is.ultrametric can be helpful. Why not simply try e.g.
>>> is.ultrametric(tree, tol=.01) 
>>> So in this sense RTFM
>>> Regards,
>>> Klaus
>>>
>>> On Aug 16, 2016 9:31 AM, "Martin Dohrmann" <
>>> m.dohrm...@lrz.uni-muenchen.de>
>>> wrote:
>>>
>>>
>>>> Am 16.08.2016 um 15:20 schrieb Joseph W. Brown:
>>>>
>>>> I agree that it is almost certainly numerical precision, as rescaling
>>>>>
>>>> the tree fixes things:
>>>>
>>>>>
>>>>> is.ultrametric(phy);
>>>>>>
>>>>> [1] FALSE
>>>>>
>>>>>> fie <- phy;
>>>>>> fie$edge.length <- fie$edge.length * 0.1;
>>>>>> is.ultrametric(fie);
>>>>>>
>>>>> [1] TRUE
>>>>>
>>>>> I’ve also tested the ultrametricity(?) with a non-R program and the
>>>>>
>>>> results are the same. So can we assume this is not a PhyloBayes issue,
>>>> an

Re: [R-sig-phylo] Ultrametric tree not recognized

[Klaus Schliep]

We should consider to improve is.ultrametric to use as tolerance a certain
number significant figures (e.g. 8-10 digits). This way tol would only
depend on the tip-to-root distance and not a machine dependent minimal
value. I find the use of significant figures suits really well as a
stopping criterion for example in my ML functions, especially as it works
for small and large trees equally well.

write.tree has an option to define the number of digits, it will be
generally less than the internal representation and I suspect one could
possible export with higher precision from other programs, but will result
in larger file sizes which one has to put in consideration for MCMC /
bootstrap tree files.
As a side note on the developer site of the r-project.org homepage (not
CRAN) is a link to this nice ancient paper (
http://www.validlab.com/goldberg/paper.pdf)

Regards,
Klaus


PS: And Martin of course I am happy with Liam's solution as it uses
phangorn ;)

On Tue, Aug 16, 2016 at 10:22 AM, Joseph W. Brown <josep...@umich.edu>
wrote:

> That’s a good tip, but the issue here is that many packages use
> is.ultrametric internally with a hard-coded (likely, default) value of tol.
> In this case, Martin simply cannot perform the analyses he wants to run
> because this option is inaccessible.
>
> 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 16 Aug, 2016, at 10:15, Klaus Schliep <klaus.schl...@gmail.com> wrote:
>
> Hello all,
> this may come be surprising to many, but consulting the manual
> ?is.ultrametric can be helpful. Why not simply try e.g.
> is.ultrametric(tree, tol=.01) 
> So in this sense RTFM
> Regards,
> Klaus
>
> On Aug 16, 2016 9:31 AM, "Martin Dohrmann" <m.dohrm...@lrz.uni-muenchen.de>
> wrote:
>
>>
>> Am 16.08.2016 um 15:20 schrieb Joseph W. Brown:
>>
>> > I agree that it is almost certainly numerical precision, as rescaling
>> the tree fixes things:
>> >
>> > > is.ultrametric(phy);
>> > [1] FALSE
>> > > fie <- phy;
>> > > fie$edge.length <- fie$edge.length * 0.1;
>> > > is.ultrametric(fie);
>> > [1] TRUE
>> >
>> > I’ve also tested the ultrametricity(?) with a non-R program and the
>> results are the same. So can we assume this is not a PhyloBayes issue, and
>> rather an unavoidable problem with large, old trees? But the fact that
>> MrBayes trees are ultrametric is confusing; does MrBayes use less precise
>> edge lengths?
>>
>> Actually the MrBayes branch lengths appear to be MORE precise.
>>
>> Martin
>>
>> >>
>> >> On 16 Aug, 2016, at 08:53, Liam J. Revell <liam.rev...@umb.edu> wrote:
>> >>
>> >> Hi Martin.
>> >>
>> >> Since you are writing & reading trees to file, my guess is that it has
>> to do with numerical precision - that is, the rounding of your edge lengths
>> when they are written to file.
>> >>
>> >> Does your tree look ultrametric when plotted in R? If so, this is
>> probably the case.
>> >>
>> >> My recommendation is that you use phangorn to compute the non-negative
>> least-squares edge lengths with the condition that the tree is ultrametric.
>> This will give you the edge lengths that result in the distances between
>> taxa with minimum sum of squared differences from the distances implied by
>> your input tree, under the criterion that the resulting tree is ultrametric.
>> >>
>> >> To do this you need to merely run:
>> >>
>> >> library(phytools)
>> >> library(phangorn)
>> >> is.ultrametric(tree) ## fails
>> >> plotTree(tree,ftype="off") ## does my tree look ultrametric?
>> >> nnls<-nnls.tree(cophenetic(tree),tree,rooted=TRUE)
>> >> is.ultrametric(tree) ## should pass
>> >>
>> >> Let us know if this works. All the best, Liam
>> >>
>> >> Liam J. Revell, Associate Professor of Biology
>> >> University of Massachusetts Boston
>> >> web: http://faculty.umb.edu/liam.revell/
>> >> email: liam.rev...@umb.edu
>> >> blog: http://blog.phytools.org
>> >>
>> >> On 8/16/2016 6:41 AM, Martin Dohrmann wrote:
>> >>> Hi,
>> >>> I want to do some diversification pattern analyses with various R
>> packages. I'm using a 

Re: [R-sig-phylo] Ultrametric tree not recognized

[Klaus Schliep]

Hello all,
this may come be surprising to many, but consulting the manual
?is.ultrametric can be helpful. Why not simply try e.g.
is.ultrametric(tree, tol=.01) 
So in this sense RTFM
Regards,
Klaus

On Aug 16, 2016 9:31 AM, "Martin Dohrmann" 
wrote:

>
> Am 16.08.2016 um 15:20 schrieb Joseph W. Brown:
>
> > I agree that it is almost certainly numerical precision, as rescaling
> the tree fixes things:
> >
> > > is.ultrametric(phy);
> > [1] FALSE
> > > fie <- phy;
> > > fie$edge.length <- fie$edge.length * 0.1;
> > > is.ultrametric(fie);
> > [1] TRUE
> >
> > I’ve also tested the ultrametricity(?) with a non-R program and the
> results are the same. So can we assume this is not a PhyloBayes issue, and
> rather an unavoidable problem with large, old trees? But the fact that
> MrBayes trees are ultrametric is confusing; does MrBayes use less precise
> edge lengths?
>
> Actually the MrBayes branch lengths appear to be MORE precise.
>
> Martin
>
> >>
> >> On 16 Aug, 2016, at 08:53, Liam J. Revell  wrote:
> >>
> >> Hi Martin.
> >>
> >> Since you are writing & reading trees to file, my guess is that it has
> to do with numerical precision - that is, the rounding of your edge lengths
> when they are written to file.
> >>
> >> Does your tree look ultrametric when plotted in R? If so, this is
> probably the case.
> >>
> >> My recommendation is that you use phangorn to compute the non-negative
> least-squares edge lengths with the condition that the tree is ultrametric.
> This will give you the edge lengths that result in the distances between
> taxa with minimum sum of squared differences from the distances implied by
> your input tree, under the criterion that the resulting tree is ultrametric.
> >>
> >> To do this you need to merely run:
> >>
> >> library(phytools)
> >> library(phangorn)
> >> is.ultrametric(tree) ## fails
> >> plotTree(tree,ftype="off") ## does my tree look ultrametric?
> >> nnls<-nnls.tree(cophenetic(tree),tree,rooted=TRUE)
> >> is.ultrametric(tree) ## should pass
> >>
> >> Let us know if this works. All the best, Liam
> >>
> >> Liam J. Revell, Associate Professor of Biology
> >> University of Massachusetts Boston
> >> web: http://faculty.umb.edu/liam.revell/
> >> email: liam.rev...@umb.edu
> >> blog: http://blog.phytools.org
> >>
> >> On 8/16/2016 6:41 AM, Martin Dohrmann wrote:
> >>> Hi,
> >>> I want to do some diversification pattern analyses with various R
> packages. I'm using a time-calibrated tree produced by PhyloBayes. Branch
> lengths are in millions of years and all taxa are extant, so the tree
> should be ultrametric. However, when I call "is.ultrametric", it returns
> "FALSE".
> >>>
> >>> Has anybody encountered something like this? Any ideas about what's
> going on/how to solve this?
> >>>
> >>> Some further information:
> >>> I also tried other PhyloBayes time trees, with the same result. In
> contrast, MrBayes time trees I tried for comparison are recognized as
> ultrametric. Regarding my diversification analyses, TESS would not run,
> telling me "The likelihood function is only defined for ultrametric
> trees!". On the other hand, BAMMtools doesn't seem to have a problem with
> my tree. I haven't tried other packages yet, but I suspect RPANDA, TreePar
> etc. might also have issues if they don't recognize my tree as ultrametric.
> >>>
> >>> I'd appreciate any help!
> >>>
> >>> Best wishes,
> >>> Martin
> >>>
> >>> Dr. Martin Dohrmann
> >>> Ludwig-Maximilians-University Munich
> >>> Dept. of Earth & Environmental Sciences
> >>> Palaeontology & Geobiology
> >>> Molecular Geo- & Palaeobiology Lab
> >>> Richard-Wagner-Str. 10
> >>> 80333 Munich, Germany
> >>> Phone: +49-(0)89-2180-6593
> >>>
> >>>
> >>> [[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-ph...@r-project.org/
> >>>
> >>
> >> ___
> >> 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-ph...@r-project.org/
> >
>
> Dr. Martin Dohrmann
> Ludwig-Maximilians-University Munich
> Dept. of Earth & Environmental Sciences
> Palaeontology & Geobiology
> Molecular Geo- & Palaeobiology Lab
> Richard-Wagner-Str. 10
> 80333 Munich, Germany
> Phone: +49-(0)89-2180-6593
>
>
> [[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-ph...@r-project.org/
>

[[alternative HTML version deleted]]

___
R-sig-phylo mailing list - 

[R-sig-phylo] Nantucket phylogeny developeR bootcamp

[Klaus Schliep]

Nantucket phylogeny developeR bootcamp

We are happy to announce a new graduate-level workshop on phylogenetic
method development R. The course will be five days in length and will
take place at the University of Massachusetts Boston's Nantucket Field
Station from the 22nd to the 26th of August, 2016. This workshop is
primarily intended for evolutionary biologists with some prior
experience in computer programming who have an interest in participating
more fully in the community of method developers for the R statistical
computing environment in phylogenetics.

On the first 1.5 days of the workshop, course leaders Drs. April Wright
and Klaus Schliep (with Dr. Liam Revell participating remotely) will
provide an introduction to the primary data structures and methods of
common phylogenetic R packages, the basics of computational algorithms
for phylogenies, and an overview of package development in R, along with
other essential topics that will depend on the prior experience of the
bootcamp participants. Over subsequent workshop days participants will
work in groups to develop small R packages on their chosen topics. These
projects will focus on adding new functionality to existing R software
in phylogenetics, and might range from tree manipulation, to
phylogenetic inference, to comparative methods, to phylogenetic
simulations, to the visualization of phylogenies or macroevolutionary
data on trees.

The workshop is funded by a grant from the National Science Foundation
to Dr. Liam Revell, with additional support from the University of
Massachusetts Boston. All accepted students will be offered a stipend to
cover or defray travel costs to and from Nantucket, and room and board
during the workshop will be provided. As the workshop will be held at a
field station, accommodation is comfortable, but basic, and participants
should be prepared to stay in multiple occupancy rooms.

To apply for the course, please fill the google form
(http://goo.gl/forms/4m5ILAGfzRihe9jS2) and submit your CV along with a
short (1 page) description of your research interests and your reasons
for participating in the workshop. Details of relevant programming
background (computer languages, R knowledge, GitHub repositories, ...)
should also be provided. Admission is competitive, and preference will
go towards students with background in evolutionary biology, basic to
moderate experience in computer programming, and a compelling motivation
for taking the course. Applications should be submitted the google form
(http://goo.gl/forms/4m5ILAGfzRihe9jS2) by July 15th, 2016. Questions
can be directed to klaus.schl...@umb.edu.

For an HTML version of this advertisement, please see:
http://www.phangorn.org/pages/Nantucket.html



-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston

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

Re: [R-sig-phylo] Error in plotBS when number of taxa is small

[Klaus Schliep]

HI Kamilla,
can you send me your code and data and what phangorn version are you using?
Cheers,
Klaus

On Tue, Jun 14, 2016 at 6:42 PM, Kamila Naxerova 
wrote:

> Hi all,
>
> I am trying to attach bootstrap values to a small tree (just 4 taxa).
> plotBS() in the phangorn package gives me the following error:
>
> tree <- plotBS(phylotree,bstrees,type="phylogram",p=90)
>
> Error in if (drop[j]) next : missing value where TRUE/FALSE needed
>
>
> This error disappears and everything work great as soon as I have 5 taxa.
> But removing one always leads to this error, for multiple independent trees.
>
> Am I missing something obvious? Unfortunately I don’t understand the error
> message.
>
> Many thanks!
> Kamila
>
>
>
> The information in this e-mail is intended only for th...{{dropped:28}}

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

Re: [R-sig-phylo] new testing version of ape

[Klaus Schliep]

Hi Emmanuel,
agree looks cool. When do you want to push the new version to CRAN? Newest
phangorn development has now a as.AAbin.phyDat function and dist.ml has
transition rate matrices for amino acids.
Does anybody know how to set .travis.yaml on github to install ape from
Emmanuel's repository instead of CRAN server?
Cheers,
Klaus

On Mon, Feb 22, 2016 at 5:37 AM, Thibaut Jombart <thibautjomb...@gmail.com>
wrote:

> Hi Emmanuel,
>
> this looks really cool - especially the new AAbin class. Looking forward to
> playing with the new version!
>
> Best
> Thibaut
>
>
> --
> Dr Thibaut Jombart
> Lecturer, Department of Infectious Disease Epidemiology
> Imperial College London
> https://sites.google.com/site/thibautjombart/
> https://github.com/thibautjombart
> Twitter: @TeebzR <https://twitter.com/TeebzR>
>
> On 22 February 2016 at 09:16, Emmanuel Paradis <emmanuel.para...@ird.fr>
> wrote:
>
> > Dear all,
> >
> > A new testing version of ape (3.4-0.3) is available. It includes three
> > main new features:
> >
> > - A new data class, "AAbin", to store amino acid sequences; there are
> > eleven new functions to generate and manipulate them including
> translation
> > from DNA.
> >
> > - The function checkAlignment does some diagnostics on a DNA alignment.
> >
> > - Some extra plotting functions built on ape's basic plotting functions
> > (plot.phylo, nodelabels, phydataplot, ring) for a specific task. For the
> > moment, there are two such functions: drawSupportOnEdges() and
> > plotBreakLongEdges().
> >
> > Information on the new functions can be found at:
> >
> > http://ape-package.ird.fr/NEWS
> >
> > And instructions on how to get and install the testing version here:
> >
> > http://ape-package.ird.fr/ape_installation.html#versions
> >
> > A Windows version is available.
> >
> > Best,
> >
> > Emmanuel
> >
> > ___
> > 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/
>



-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston

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

Re: [R-sig-phylo] colouring tree branches like tips

[Klaus Schliep]

Dear Gudrun,

almost right, try this code:
set.seed(42)
tree<-pbtree(n=10)
tree$tip.label<-c("A","B","Q","M","C","D","F","X","Z","W")
color<-rep("black",length(tree$tip.label))
color[which(tree$tip.label=="A")]<-"red"
edgecolor<-rep("black",length(tree$edge))
plot(tree)
edgelabels()
# edge 3 leads to A
edgecolor[3]<-"red"
plot(tree, edge.color=edgecolor)

Regards,
Klaus


On Tue, Jan 26, 2016 at 3:06 AM, Gygli, Gudrun <gudrun.gy...@wur.nl> wrote:

> Dear fellow phylo-R users,
>
>
> I am trying to colour tree (terminal and non terminal) branches in the
> same colour as the tip of said branch.
>
>
> To illustrate, see the code below, which is almost doing what I want. You
> can see that only colouring a specific tip, which is called A in this
> example, is no problem.
>
>
> However, I cannot colour the branches leading up to tip "A" in red with
> this code. I tried using tree$edge, but this is not working the way I
> thought.
>
>
> Does anyone know how to colour all the branches leading up to tip "A" in
> red?
>
>
> tree<-pbtree(n=10)
> tree$tip.label<-c("A","B","Q","M","C","D","F","X","Z","W")
> color<-rep("black",length(tree$tip.label))
> color[which(tree$tip.label=="A")]<-"red"
> edgecolor<-rep("black",length(tree$edge))
> edgecolor[which(tree$tip.label=="A")]<-"red"
>
>
> Thanks a lot for any help.
>
>
> Best regards
>
>
> G
>
>
>
> Gudrun Gygli, MSc
>
> PhD candidate
>
> Wageningen University
> Laboratory of Biochemistry
> Dreijenlaan 3
> 6703 HA Wageningen
> The Netherlands
>
> Phone  31 317483387
> e-mail: gudrun.gy...@wur.nl
>
> - - - - - - - - - - - - - - - - - -
>
> Project information:
> http://www.wageningenur.nl/en/show/Bioinformatics-structural-biology-and-molecular-modeling-of-Vanillyl-Alcohol-Oxidases-VAOs.htm
>
> ___
> 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/
>



-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston

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

Re: [R-sig-phylo] colouring tree branches like tips

[Klaus Schliep]

So this should work:

library(phytools)
set.seed(42)
tree<-pbtree(n=10)
tree$tip.label<-c("A","B","Q","M","C","D","F","X","Z","W")
color<-rep("black",length(tree$tip.label))
color[which(tree$tip.label=="A")]<-"red"
edgecolor<-rep("black",length(tree$edge))
ind = which(tree$tip.label=="A") #
edgecolor[match(ind, tree$edge[,2])]<-"red"
plot(tree, edge.color=edgecolor)

Regards,
Klaus

On Tue, Jan 26, 2016 at 7:03 AM, Gygli, Gudrun <gudrun.gy...@wur.nl> wrote:

> Dear Klaus,
>
> Dear phylo-R users,
>
>
> your code does work but is still not doing exactly what I want.
>
>
> I am hoping to get the number of the edge which leads to A automatically
> through sth similar to what I have for the colouring of the tips.
>
> (color[which(tree$tip.label=="A")]<-"red")
>
>
>
> Also, I would like to colour the other branches (from earlier "split-off-
> nodes") leading to A in red too.
>
>
> Let me try to clarify:
>
> I work on creating a tree for a protein family and want to colour branches
> based on a selection of protein names
>
> (by creating a vector which contains all the names of proteins which share
> properties),
>
> (for example using sth like this
> color[which(tree$tip.label%in%VECTORWITHNAMES)]<-"red")
>
>  and these proteins cluster in the same part of the tree, thus I want to
> colour all the branches leading to these protein-tips red (or any other
> colour).
>
> To add some complexity, there are different properties that different
> proteins in that tree share, so in the end I would like to have different
> colours in the tree, with each colour reflecting another property...
>
>
> Any ideas if this is possible in R or with another tool?
>
>
> Thanks for the help
>
>
> G
>
>
> Gudrun Gygli, MSc
>
> PhD candidate
>
> Wageningen University
> Laboratory of Biochemistry
> Dreijenlaan 3
> 6703 HA Wageningen
> The Netherlands
>
> Phone  31 317483387
> e-mail: gudrun.gy...@wur.nl
>
> - - - - - - - - - - - - - - - - - -
>
> Project information:
> http://www.wageningenur.nl/en/show/Bioinformatics-structural-biology-and-molecular-modeling-of-Vanillyl-Alcohol-Oxidases-VAOs.htm
> 
> From: Klaus Schliep <klaus.schl...@gmail.com>
> Sent: Tuesday, January 26, 2016 10:34 AM
> To: Gygli, Gudrun
> Cc: r-sig-phylo@r-project.org
> Subject: Re: [R-sig-phylo] colouring tree branches like tips
>
> Dear Gudrun,
>
> almost right, try this code:
> set.seed(42)
> tree<-pbtree(n=10)
> tree$tip.label<-c("A","B","Q","M","C","D","F","X","Z","W")
> color<-rep("black",length(tree$tip.label))
> color[which(tree$tip.label=="A")]<-"red"
> edgecolor<-rep("black",length(tree$edge))
> plot(tree)
> edgelabels()
> # edge 3 leads to A
> edgecolor[3]<-"red"
> plot(tree, edge.color=edgecolor)
>
> Regards,
> Klaus
>
>
> On Tue, Jan 26, 2016 at 3:06 AM, Gygli, Gudrun <gudrun.gy...@wur.nl
> <mailto:gudrun.gy...@wur.nl>> wrote:
> Dear fellow phylo-R users,
>
>
> I am trying to colour tree (terminal and non terminal) branches in the
> same colour as the tip of said branch.
>
>
> To illustrate, see the code below, which is almost doing what I want. You
> can see that only colouring a specific tip, which is called A in this
> example, is no problem.
>
>
> However, I cannot colour the branches leading up to tip "A" in red with
> this code. I tried using tree$edge, but this is not working the way I
> thought.
>
>
> Does anyone know how to colour all the branches leading up to tip "A" in
> red?
>
>
> tree<-pbtree(n=10)
> tree$tip.label<-c("A","B","Q","M","C","D","F","X","Z","W")
> color<-rep("black",length(tree$tip.label))
> color[which(tree$tip.label=="A")]<-"red"
> edgecolor<-rep("black",length(tree$edge))
> edgecolor[which(tree$tip.label=="A")]<-"red"
>
>
> Thanks a lot for any help.
>
>
> Best regards
>
>
> G
>
>
>
> Gudrun Gygli, MSc
>
> PhD candidate
>
> Wageningen University
> Laboratory of Biochemistry
> Dreijenlaan 3
> 6703 HA Wageningen
> The Netherlands
>
> Phone  31 317483387
> e-mail: gudrun.gy...@wur.nl<mailto:gudrun.gy...@wur.nl>
>
> - - - - - - - - - - - - - - - - - -
>
> Project information:
> http://www.wageningenur.nl/en/show/Bioinformatics-structural-biology-and-molecular-modeling-of-Vanillyl-Alcohol-Oxidases-VAOs.htm
>
> ___
> R-sig-phylo mailing list - R-sig-phylo@r-project.org 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/
>
>
>
> --
> Klaus Schliep
> Postdoctoral Fellow
> Revell Lab, University of Massachusetts Boston
>
>


-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston

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

Re: [R-sig-phylo] extract parent-child from newick

[Klaus Schliep]

Dear Lev,

I hope this little example explains it:

tree <- read.tree(text = "((t3:1,t2:1)b:1,t1:2)a;")
par(mfrow = c(1,2))
plot(tree, show.node.label=TRUE)
plot(tree, show.tip.label=FALSE)
nodelabels()
tiplabels()
tree$tip.label
tree$node.label

The tree$tip.label correspond to 1 to the number of tip labels and
tree$node.label to (number of tip labels +1) to (number of nodes) in the
tree$edge matrix.

Regards,
Klaus

On Tue, Jan 12, 2016 at 7:54 PM, Yampolsky, Lev <yampo...@mail.etsu.edu>
wrote:

> A very simple question: I have a tree with all tips and nodes labeled and
> I want parent-child pairs from it. Looking at tree$edge does not quite
> help, because these parent and child pairs are labeled by the internal
> index, not IDs from the newick file. Looking at tree$tip.label and
> tree$node.label does not help, because these two vectors index tips and
> nodes separately, so these indexes do not correspond to those in tree$edge.
>
> Should be very very obvious how to do it!
>
> Thanks!
>
> --
> Lev Yampolsky
>
> Professor
> Department of Biological Sciences
> East Tennessee State University
> Box 70703
> Johnson City TN 37614-1710
> Cell 423-676-7489
> Office/lab 423-439-4359
> Fax423-439-5958
>
> ___
> 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/
>



-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston

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

Re: [R-sig-phylo] Package phytools load failed

[Klaus Schliep]

Hi all,
you have to install the Biostrings package from bioconductor. Than it
should works. Phangorn has a dependency on biostrings and phytools depend
on phangorn.
Have a nice weekend
Klaus
Am 09.01.2016 16:35 schrieb "Donald Miles" :

> Dear All,
>
> I receive the same error as Gabriela when I try to load phytools in R or R
> Studio.
>
> Here is the text:
>
> library("phytools",
> lib.loc="/Library/Frameworks/R.framework/Versions/3.2/Resources/library")
> Error in loadNamespace(j <- i[[1L]], c(lib.loc, .libPaths()), versionCheck
> = vI[[j]]) :
>   there is no package called ‘Biostrings’
> Error: package or namespace load failed for ‘phytools’
>
> I am running R v 3.2.3 and phytools v 0.5-10.
>
> If I try to load Biostrings I receive an error message that Biostrings is
> not compatible with R 3.2.3.
>
> Cheers,
>
> Don Miles
> Ohio University
>
>
>
> On Sat, Jan 9, 2016 at 12:50 PM, Vojtěch Zeisek  wrote:
>
> > Hello
> >
> > Dne So 9. ledna 2016 13:11:45, Gabriela Wofkova napsal(a):
> > > ​Hello,
> > >
> > > I have an problem with loading package "phytools". I have the newest
> > > version of R studio​ and it cannot load package "Biostrings" which is
> > > reqired to package "phytools". Also it not works in R program, too. Can
> > you
> > > have any idea, how to evade this?
> > >
> > > I need to use "read.newick" and load phylogenetic tree.
> >
> > What exactly are You doing? phytools does not require Biostrings and
> > Biostrings does not require phytools. :-)
> > Which operating system are You running? Which error did You get? Errors
> are
> > keys to solutions... Do You have newest versions of the packages (I'd try
> > to
> > reinstall them first to see if they weren't corrupted)?
> > Both load fine for me (on Linux):
> >
> > > library(phytools)
> > Loading required package: ape
> > Loading required package: maps
> >
> >  # ATTENTION: maps v3.0 has an updated 'world' map.#
> >  # Many country borders and names have changed since 1990. #
> >  # Type '?world' or 'news(package="maps")'. See README_v3. #
> >
> > > library(Biostrings)
> > Loading required package: BiocGenerics
> > Loading required package: parallel
> >
> > Attaching package: ‘BiocGenerics’
> >
> > The following objects are masked from ‘package:parallel’:
> >
> > clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
> clusterExport,
> > clusterMap, parApply, parCapply, parLapply, parLapplyLB,
> > parRapply, parSapply, parSapplyLB
> >
> > The following objects are masked from ‘package:stats’:
> >
> > IQR, mad, xtabs
> >
> > The following objects are masked from ‘package:base’:
> >
> > anyDuplicated, append, as.data.frame, as.vector, cbind, colnames,
> > do.call,
> > duplicated, eval, evalq, Filter, Find, get, grep, grepl,
> > intersect, is.unsorted, lapply, lengths, Map, mapply, match, mget,
> > order,
> > paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
> > rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply, union,
> > unique, unlist, unsplit
> >
> > Loading required package: S4Vectors
> > Loading required package: stats4
> > Loading required package: IRanges
> > Loading required package: XVector
> >
> > You do not need Biostrings to use read.newick.
> > Note there is also read.tree in ape package reading trees in Newick
> format.
> >
> > > ​Thanks a lot for any useful answer!​
> > >
> > > Gabriela, Charles University in Prague​
> >
> > Sincerely,
> > Vojtěch
> >
> > --
> > Vojtěch Zeisek
> > http://trapa.cz/en/
> >
> > Department of Botany, Faculty of Science
> > Charles University in Prague
> > Benátská 2, Prague, 12801, CZ
> > http://botany.natur.cuni.cz/en/
> >
> > Institute of Botany, Academy of Science
> > Zámek 1, Průhonice, 25243, CZ
> > http://www.ibot.cas.cz/en/
> >
> > Czech Republic
> >
> > ___
> > 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/

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

Re: [R-sig-phylo] Calculating branch lengths on a multiphylo object

[Klaus Schliep]

Happy New Year,
it seems you have not properly installed the packages and you may want also
update R.
Cheers, Klaus

On Fri, Jan 1, 2016 at 12:34 AM, Michael Pauers <mjpau...@gmail.com> wrote:

> Dear list,
>
> First, Happy New Year!
>
> Second, I will apologize in advance, as I am nothing more than a duffer
> with R, so I may ask several follow-up questions to this one.
>
> I have a set of 1000 trees that lack branch lengths.  I have used
> previously used ape to calculate branch lengths on a single tree, but never
> on a multiphylo object.  I did find and have tried the approach documented
> at this link:
>
> http://blog.phytools.org/2015/04/phylogenetic-regression-when-branch.html
>
> but I received a few error messages, including:
>
> > library(phytools)
> Error in gzfile(file, "rb") : cannot open the connection
> In addition: Warning messages:
> 1: package ‘phytools’ was built under R version 3.0.3
> 2: In read.dcf(file.path(p, "DESCRIPTION"), c("Package", "Version")) :
>   cannot open compressed file 'C:/Program
> Files/R/R-3.0.2/library/ape/DESCRIPTION', probable reason 'No such file or
> directory'
> 3: In gzfile(file, "rb") :
>   cannot open compressed file '', probable reason 'No such file or
> directory'
> > library(phangorn)
> Error in gzfile(file, "rb") : cannot open the connection
> In addition: Warning messages:
> 1: package ‘phangorn’ was built under R version 3.0.3
> 2: In read.dcf(file.path(p, "DESCRIPTION"), c("Package", "Version")) :
>   cannot open compressed file 'C:/Program
> Files/R/R-3.0.2/library/ape/DESCRIPTION', probable reason 'No such file or
> directory'
> 3: In gzfile(file, "rb") :
>   cannot open compressed file '', probable reason 'No such file or
> directory'
>
> and:
>
> > library(ape)
> > ## simulate 200 pure-birth trees:
> > trees<-ttrees<-pbtree(n=50,scale=1,nsim=200)
> Error: could not find function "pbtree"
> > ## ttrees contains the trees with their original branch lengths
> > ## we'll perform manipulations on trees
> > foo<-function(tree){
> + obj<-fastBM(tree,nsim=2)
> + colnames(obj)<-c("x","y")
> + as.data.frame(obj)
> + }
> > xy<-lapply(ttrees,foo)
> Error in lapply(ttrees, foo) : object 'ttrees' not found
> > library(nlme)
> > fit.model<-function(tree,data){
> + data$v<-diag(vcv.phylo(tree))
> + fit<-gls(y~x,data=data,correlation=corBrownian(1,tree),
> + weights=varFixed(~v))
> + setNames(c(coefficients(fit)[2],anova(fit)$"p-value"[2]),
> + c("beta","p-value"))
> + }
> > fit.true<-t(mapply(fit.model,trees,xy))
> Error in mapply(fit.model, trees, xy) : object 'xy' not found
> > mean(fit.true[,1]) ## should be zero
> Error in mean(fit.true[, 1]) : object 'fit.true' not found
> > noogtree<-lapply(noogtree,compute.brlen)
> > class(noogtree)<-"multiPhylo"
> > fit.grafen<-t(mapply(fit.model,noogtree,xy))
> Error in mapply(fit.model, noogtree, xy) : object 'xy' not found
> > utils:::menuInstallPkgs()
> Error in install.packages(NULL, .libPaths()[1L], dependencies = NA, type =
> type) :
>   no packages were specified
> > fit.grafen<-t(mapply(fit.model,noogtree))
> Error in data$v <- diag(vcv.phylo(tree)) :
>   argument "data" is missing, with no default
> > fit.grafen<-t(mapply(fit.model,noogtree,xy))
> Error in mapply(fit.model, noogtree, xy) : object 'xy' not found
> > mean(fit.grafen[,1])
> Error in mean(fit.grafen[, 1]) : object 'fit.grafen' not found
> > foo<-function(tree){
> + tree$edge.length<-rep(1,nrow(tree$edge))
> + tree
> + }
> > trees<-lapply(ttrees,foo)
> Error in lapply(ttrees, foo) : object 'ttrees' not found
> > class(trees)<-"multiPhylo"
> > fit.equal<-t(mapply(fit.model,trees,xy))
> Error in mapply(fit.model, trees, xy) : object 'xy' not found
> > mean(fit.equal[,1]) ## should be zero
> Error in mean(fit.equal[, 1]) : object 'fit.equal' not found
>
> Does anyone have any advice?  Please let me know if I should provide more
> information.
>
> Thanks in advance,
>
> Mike Pauers
>
> [[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/




-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston

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

Re: [R-sig-phylo] Detecting homoplasy

[Klaus Schliep]

Dear Tim,

I added an optional argument sitewise to the CI and RI function in
phangorn. The vector may contains NaN if the function is undefined: pscore
of 0 for CI, uninformative site for RI.
It is now on github:
library(devtools)
install_github("KlausVigo/phangorn")

Cheers,
Klaus


On Thu, Dec 17, 2015 at 8:38 PM, Read, Timothy D  wrote:

> Hi all,
>
>
> I am interested in finding homoplasious sites in a nucleotide sequence
> alignment on a ML tree calculated from that alignment.  Is there a commonly
> used approach for doing this in ape, phangorn etc.  I have tried a number
> of google searches based around 'R homoplasy'  but cant find anything.
>
>
> On a related note, is there a function available to calculate the
> consistency index for a given position on an alignment rather than an
> average for the whole alignment, which is what CI in phangorn does now?
>
>
> thanks!
>
>
> Tim
>
> 
>
> This e-mail message (including any attachments) is for...{{dropped:19}}

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

[R-sig-phylo] Problems with dependency Biostrings

[Klaus Schliep]

 Hello all,

there have been a few reports of problems with installing phytools or
phangorn and you may see something like this:

devtools::install_github("liamrevell/phytools")
Downloading GitHub repo liamrevell/phytools@master
Installing phytools
Skipping 1 packages not available: Biostrings
...
...
Error: Command failed (1)

Recently I updated phangorn (now version 2.0.1), but there were also
updates of R (3.2.3) and also Bioconductor. All this seem to have
contributed to this problem. install.packages() and install.github() do not
always get bioconductor dependencies right (e.g.
https://github.com/hadley/devtools/issues/700).

To solve the problem simply install the Biostrings package first

source("http://bioconductor.org/biocLite.R;)
biocLite("Biostrings")

and in a fresh session of R try again to install phytools or phangorn.

If you have installed any bioconductor packages, chances are that you also
installed Biostrings package are pretty high.

Regards,
Klaus



-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston

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

Re: [R-sig-phylo] Mapping gene gain loss info to rooted species tree

[Klaus Schliep]

Hi Annand,

I added a short script and the output of it. This should help you to get
started. Have a look in the documentation of plot.phylo to make it nicer
looking.

Cheers,
Klaus

On Mon, Nov 16, 2015 at 6:55 PM, Anand K S Rao <aks...@ucdavis.edu> wrote:

> Dear folks,
>
> I seek your help with mapping gene gain and loss numbers for leaf nodes and
> internal nodes to a rooted species tree.
>
> I have 2 pieces of information.
>
> *1. *A species tree in NEWICK format that has internal node numbering info
> apart from leaf node names:
>
> (((Mpu229:1.0,Mpu228:1.0)n2:1.0,(((Vca199:1.0,Cre169:1.0)n5:1.0,Csu227:1.0)n7:1.0,((Car:1.0,'Mt3.5v5':1.0)n10:1.0,((Gma109:1.0,Pvu186:1.0)n13:1.0,Cca:1.0)n15:1.0)n16:1.0,(((Mdo196:1.0,Ppe139:1.0)n19:1.0,Fve226:1.0)n21:1.0,Csa122:1.0)n23:1.0)n24:1.0,Aly107:1.0,Ath167:1.0)n27:1.0,Cru183:1.0)n29:1.0,(Tha173:1.0,Bra197:1.0)n32:1.0)n33:1.0,Cpa113:1.0)n35:1.0,(Tca233:1.0,Gra221:1.0)n38:1.0)n39:1.0,(Ccl165:1.0,Csi154:1.0)n42:1.0)n43:1.0,((Rco119:1.0,Mes147:1.0)n46:1.0,(Ptr156:1.0,Lus200:1.0)n49:1.0)n50:1.0)n51:1.0,Egr201:1.0)n53:1.0)n54:1.0,Vvi145:1.0)n56:1.0,((Sly225:1.0,Stu206:1.0)n59:1.0,Mgu140:1.0)n61:1.0)n62:1.0,Aco195:1.0)n64:1.0,(((Zma181:1.0,Sbi79:1.0)n67:1.0,(Pvi202:1.0,Sit164:1.0)n70:1.0)n71:1.0,(Bdi192:1.0,Osa193:1.0)n74:1.0)n75:1.0)n76:1.0,Smo91:1.0)n78:1.0,Ppa152:1.0)n80:2.0)n81:1.0)n82:0.5,Olu231:0.5)n84;
>
> and
>
> *2.* Gene gain and loss information of each of the leaf nodes and internal
> nodes (1st row - names / numbers of nodes corresponding to NEWICK tree
> above, 2nd row - Gene gain numbers, 3rd row - Gene loss numbers)
> Mpu229 n2 Mpu228 n82 Vca199 n5 Cre169 n7 Csu227 n81 Car n10 Mt3.5v5 n16
> Gma109 n13 Pvu186 n15 Cca n24 Mdo196 n19 Ppe139 n21 Fve226 n23 Csa122 n54
> Aly107 n27 Ath167 n29 Cru183 n33 Tha173 n32 Bra197 n35 Cpa113 n39 Tca233
> n38
> Gra221 n43 Ccl165 n42 Csi154 n51 Rco119 n46 Mes147 n50 Ptr156 n49 Lus200
> n53
> Egr201 n56 Vvi145 n62 Sly225 n59 Stu206 n61 Mgu140 n64 Aco195 n76 Zma181
> n67
> Sbi79 n71 Pvi202 n70 Sit164 n75 Bdi192 n74 Osa193 n78 Smo91 n80 Ppa152 n84
> Olu231
> 0 0 0 0 0 3 2 0 0 12 3 1 220 36 38 1 9 26 11 8 76 14 27 0 0 0 3 392 50 18 9
> 0 0 57 33 5 29 0 3 0 0 0 0 0 3 38 2 9 1 3 6 13 35 1 80 0 13 13 5 87 0 0 0 0
> 90 113 31 131 17 9 11 82 136 34 42 219 51 7 60 59 17 54 14 0 0
> 0 0 0 0 13 0 8 2 23 0 109 33 103 121 30 11 49 191 54 428 29 0 54 34 112 211
> 103 17 52 96 92 124 179 38 66 227 83 0 342 81 0 380 0 124 29 374 35 37 20
> 100 19 435 65 32 60 163 565 22 355 21 0 68 0 329 0 67 254 12 81 231 42 94
> 111 59 152 132 104 138 118 1 107 7 55 0 0
>
> It is possible to use any R/BioConductoR phylogeny packages (like "ape") to
> draw the species tree with the gene gains and losses mapped to the nodes?
>
> I am somewhere between an R noob and novice. So example syntax would be
> most helpful and appreciated.
>
> I hope you wont blame me for wanting to get fancy and allow coloring of
> gains and losses
>
> Thanks a ton!
> Anand
>
> --
> Anand K.S. Rao  PhD candidate, Plant Biology <http://www-plb.ucdavis.edu/>
> with a Designated Emphasis in Biotechnology <http://deb.ucdavis.edu/>,
> Doug Cook Laboratory for Legume Genetics and Genomics
> <http://cooklab.ucdavis.edu/>, Dept. of Plant Pathology
> <http://plantpathology.ucdavis.edu/>.
> UC- Davis <http://ucdavis.edu/>,  CA - 95616 USA   |   aks...@ucdavis.edu
>  |   (530) 574-5134   |   LinkedIn <http://www.linkedin.com/in/anandksrao>
> _
>   CTTATTGTTGAACTTOAATGGTGCTAATGATCCTCGTOTCTCCTGAACGT - translate THAT!
>
> [[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/
>



-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston
Mpu229 n2 Mpu228 n82 Vca199 n5 Cre169 n7 Csu227 n81 Car n10 Mt3.5v5 n16 Gma109 
n13 Pvu186 n15 Cca n24 Mdo196 n19 Ppe139 n21 Fve226 n23 Csa122 n54 Aly107 n27 
Ath167 n29 Cru183 n33 Tha173 n32 Bra197 n35 Cpa113 n39 Tca233 n38 Gra221 n43 
Ccl165 n42 Csi154 n51 Rco119 n46 Mes147 n50 Ptr156 n49 Lus200 n53 Egr201 n56 
Vvi145 n62 Sly225 n59 Stu206 n61 Mgu140 n64 Aco195 n76 Zma181 n67 Sbi79 n71 
Pvi202 n70 Sit164 n75 Bdi192 n74 Osa193 n78 Smo91 n80 Ppa152 n84 Olu231
0 0 0 0 0 3 2 0 0 12 3 1 220 36 38 1 9 26 11 8 76 14 27 0 0 0 3 392 50 18 9 0 0 
57 33 5 29 0 3 0 0 0 0 0 3 38 2 9 1 3 6 13 35 1 80 0 13 13 5 87 0 0 0 0 90 113 
31 131 17 9 11 82 136 34 42 219 51 7 60 59 17 54 14 0 0 
0 0 0 0 13 0 8 2 23 0 109 33 103 121 30 11 49 191 54 428 29 0 54 34 112 211 103 
17 52 96 92 124 1

Re: [R-sig-phylo] output of phangorn::midpoint is not compatible with ape::ladderize?

[Klaus Schliep]

Dear Guangchuang & Emmanuel,

it seems ladderize assumes that trees are in "cladewise" order, but
midpoint returns a tree in "postorder". Emmanuel can you include a
phy <- reorder(phy)
in ladderize. That should solve the problem.

x<-rtree(10)
x2 <- reorder(x, "postorder")
plot(ladderize(x2))  # error message

y<-phangorn::midpoint(x)
y <- reorder(y)
plot(ladderize(y)) # works!!!

Cheers,
Klaus



On Sun, Oct 4, 2015 at 10:22 AM, Yu, Guangchuang <g...@connect.hku.hk>
wrote:

> Dear all,
>
> In the issue, https://github.com/GuangchuangYu/ggtree/issues/15, we found
> that midpoint() is not compatible with ladderize().
>
> x<-rtree(10)
> y<-phangorn::midpoint(x)
>
> > y
> Phylogenetic tree with 10 tips and 9 internal nodes.
> Tip labels:
> t8, t1, t4, t6, t10, t3, ...
> Rooted; includes branch lengths.> y$edge
>   [,1] [,2]
>  [1,]   153
>  [2,]   154
>  [3,]   142
>  [4,]   14   15
>  [5,]   13   14
>  [6,]   135
>  [7,]   12   13
>  [8,]   126
>  [9,]   181
> [10,]   18   12
> [11,]   198
> [12,]   199
> [13,]   17   10
> [14,]   17   19
> [15,]   167
> [16,]   16   17
> [17,]   11   16
> [18,]   11   18> ape::ladderize(y)
> Phylogenetic tree with 10 tips and 9 internal nodes.
> Tip labels:
> t8, t1, t4, t6, t10, t3, ...
> Rooted; includes branch lengths.> ape::ladderize(y)$edge
>   [,1] [,2]
>  [1,]   11   18
>  [2,]   18   12
>  [3,]   16   17
>  [4,]   17   19
>  [5,]   198
>  [6,]   199
>  [7,]   154
>  [8,]   17   10
>  [9,]   181
>
> any idea?
>
>
> Bests,
> Guangchuang
>
> --
> --~--~-~--~~~---~--~~
> Guangchuang Yu, PhD Candidate
> State Key Laboratory of Emerging Infectious Diseases
> School of Public Health
> The University of Hong Kong
> Hong Kong SAR, China
> www: http://ygc.name
> -~--~~~~--~~--~--~---
>
> [[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/
>



-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston

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

Re: [R-sig-phylo] Cannot plot above 500 superimposed trees with densiTree (Phangorn)

[Klaus Schliep]

Hello Marina,

the problem seems with the alpha parameter. It seems it does not work for
smaller than 1/500, just set it to a small value, e.g. 0.01. I will include
a check into the function for further releases. The alpha value and the
edge weight are nice parameters to play with the appearance of the plot.

Regards,
Klaus

On Sat, Aug 22, 2015 at 1:51 PM, Marina Strelin marina.streli...@gmail.com
wrote:

 Dear all,


 I want to plot superimposed trees from a posterior distribution using the
 densiTree function of the R package phangorn. I’m using the R vesion 3.1.3
 (I understand that this package was written under that version). The
 function works well when I use no more than 500 trees from the posterior
 distribution, but once I reach that limit the function densiTree runs in
 the R console (no error message pops up)  but no plot is retrieved. I
 believe my computer has no memory problems, since 1/3 of the C disc is
 still empty.


 trees-read.nexus(file=file.choose())

 consensus-read.nexus(file=file.choose())

 x-sample(trees, 500)

 densiTree(x, type = cladogram, alpha = 1/length(x), consensus =
 consensus, optim = FALSE,

 scaleX = FALSE, col = 1, width = 1, cex = 0.8)


 (THIS WORKS)


 x-sample(trees, 550)

 densiTree(x, type = cladogram, alpha = 1/length(x), consensus =
 consensus, optim = FALSE,

 scaleX = FALSE, col = 1, width = 1, cex = 0.8)


 (THIS DOES NOT WORK!)


 Thanks in advance for helping me with this!


 Marina

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




-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston

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

Re: [R-sig-phylo] Plotting error edgelabels in phylogeny type=fan

[Klaus Schliep]

Dear Julia,

this function should fix your problem and Emmanuel may can include it into
edgelables function:

edgeLabelsFan - function (text, edge, adj = c(0.5, 0.5), frame = rect,
pch = NULL,
  thermo = NULL, pie = NULL, piecol = NULL, col = black,
  bg = lightgreen, horiz = FALSE, width = NULL, height = NULL,
  date = NULL, ...)
{
lastPP - get(last_plot.phylo, envir = .PlotPhyloEnv)
if (missing(edge)) {
sel - 1:dim(lastPP$edge)[1]
subedge - lastPP$edge
}
else {
sel - edge
subedge - lastPP$edge[sel, , drop = FALSE]
}
r - sqrt(lastPP$xx^2 + lastPP$yy^2)

XX - lastPP$xx[subedge[,2]] * (r[subedge[,2]] + r[subedge[,1]]) /
(r[subedge[,2]] * 2)
YY - lastPP$yy[subedge[,2]] * (r[subedge[,2]] + r[subedge[,1]]) /
(r[subedge[,2]] * 2)
if (!is.null(date))
XX[] - max(lastPP$xx) - date
BOTHlabels(text, sel, XX, YY, adj, frame, pch, thermo, pie,
   piecol, col, bg, horiz, width, height, ...)
}




Kind regards,
Klaus





On Thu, Jul 30, 2015 at 1:10 PM, Julia Dupin julia.du...@gmail.com wrote:

 Dear all,

 I'm trying to plot a tree with some branches marked using edgelabels() but
 the function is misplacing the labels (it misses the actual branches) when
 I use it in a fan type tree.
 Here is an example

 tree-pbtree(n=10,scale=10)
 plot(tree,type=fan)
 edgelabels()

 The phylogeny I am working on is pretty big so fan is the best option for
 me to visualize it.

 Has anyone had this issue before? Is there a way to fix it? Or maybe an
 alternative to edgelabels?

 Thanks in advance!!

 Julia

 --
 Julia Dupin
 PhD candidate
 Smith lab
 Dept. Ecology and Evolutionary Biology, Ramaley Hall C127A
 University of Colorado - Boulder
 email: julia.du...@colorado.edu
 website: http://www.colorado.edu/smithlab/

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




-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston

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

Re: [R-sig-phylo] Prune very large tree

[Klaus Schliep]

Have you looked into drop.tip from ape? if your labels you are interested
in a vector, tip this should work:
smalltree - drop.tip(tree, which(!(tree$tip.label %in% tip)))


Regards,
Klaus


On Wed, Jul 15, 2015 at 12:15 PM, Sergio Ferreira Cardoso 
sff.card...@campus.fct.unl.pt wrote:

 Dear members,

 I have a tree with almost 10.000 species and I want to prune it. The
 problem is I only want a tree with 59 of those species. Is there any
 possible way to make this in R? Is there any code to prune all taxa except
 the ones on my data?

 Thanks in advance.

 Best regards,
 Sérgio.

 --
 Com os melhores cumprimentos,
 Sérgio Ferreira Cardoso.

 

 Best regards,
 Sérgio Ferreira Cardoso




 MSc. Paleontology candidate
 Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa
 Geociências - Universidade de Évora

 Lisboa, Portugal
 ᐧ

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




-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston

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

Re: [R-sig-phylo] Bug in reorder.phylo() (related to cleaning phylo objects)

[Klaus Schliep]

Hi David,

I found the bug. Somehow ape assumes on the one hand that the root is
number of tips + 1.
On the other hand the root of an phylo object is the node which is in
tree$edge[,1], but not in tree$edge[,2],
this is the case even in an unrooted tree. This annoys me for a long time.

The root of tree1 with the definition above is actually node 151 and not
131. Changing these solves your problem:
Here is a small function to do this for you:

switch.nodes- function(tree, a, b){
ntips = length(tree$tip.label)
tree$edge[tree$edge==a] = 0L
tree$edge[tree$edge==b] = as.integer(a)
tree$edge[tree$edge==0L] = as.integer(b)
if(!is.null(tree$node.label)){
tmp-tree$node.label[a-ntips]
tree$node.label[a-ntips] = tree$node.label[b-ntips]
tree$node.label[b-ntips] = tmp
}
tree
}

library(phangorn)
attr(tree1, order) = NULL
getRoot(tree1)

tree2 = switch.nodes(tree1, 131, 151)
plot(tree2) # now works for me

Cheers and have a nice weekend,
Klaus







On Fri, Jun 12, 2015 at 2:53 PM, David Bapst dwba...@gmail.com wrote:

 Hello all,

 As those of you who directly manipulate the guts of phylo objects in
 your code (or construct new phylo objects whole cloth from
 un-phylo-like data structures) have probably experienced, it is
 sometimes easy to create $edge matrices that aren't accepted by ape
 functions (I often use plot.phylo as my litmus for this).

 When this occurs, various things can be done; I usually do the
 following sequence:

 collapse.singles()
 reorder.phylo(,cladewise)
 read.tree(text=write.tree(,file=NULL))

 ...or something along those lines.

 However, I've run into an issue where reorder.phylo() returns what is
 essentially a scrap of the original $edge matrix, without warning.
 Very worrisome! I'm using R version 3.2.0 and ape 3.3. Here's, my
 script, including a call to a saved object on my website, so that the
 error can be reproduced:

 (Why do I have an .Rdata object listed with a .txt extension, you
 might wonder? Because my school apparently sanitizes file extensions
 on our hosted websites that it thinks are executables, archives or
 doesn't recognize, but doesn't bother to check file innards when it
 does recognize the extensions. Its easily hacked, at least.)

 #

 library(ape)

 load(url(http://webpages.sdsmt.edu/~dbapst/weirdTree_06-12-15.txt;))

 #can I plot it?
 plot(tree)
 #nope

 tree1-collapse.singles(tree)

 #any single-child nodes left?
 sum(sapply(unique(tree1$edge[,1]),function(x) sum(x==tree1$edge[,1])==1))
 #nope

 #can I plot it?
 plot(tree1)
 #nope

 #now reorder
 tree2-reorder.phylo(tree1,cladewise)
 tree2$edge

 #now reorder with postorder
 tree2-reorder.phylo(tree1,postorder)
 tree2$edge

 #

 As we can see, reorder.phylo with various ordering methods returns an
 edge matrix with just six rows, from a tree that originally had
 hundreds of edges. Most worrisome, it does this without any error
 message. In this particular case, I retain the tree correctly if I
 skip reorder.phylo and use the read.tree(write.tree) trick, but this
 behavior of reorder.phylo is worrying.

 Anyone have a clue what's going on here? If I am perhaps misusing
 reorder.phylo(), is there an alternative approach for use in cleaning
 phylo objects?

 Cheers,
 -Dave

 --
 David W. Bapst, PhD
 Adjunct Asst. Professor, Geology and Geol. Eng.
 South Dakota School of Mines and Technology
 501 E. St. Joseph
 Rapid City, SD 57701

 http://webpages.sdsmt.edu/~dbapst/
 http://cran.r-project.org/web/packages/paleotree/index.html

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




-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston

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

Re: [R-sig-phylo] Phylogenetic analysis with sequencing error

[Klaus Schliep]

Dear George

after all this useful background information from Joe, Brian et al. back to
your initial question, if it is still valid by now. Yes, you can use pml
for this. You actually do not have to mess with the pml function itself.

As nucleotides, amino acids, codons etc. are discrete phyDat objects in
phangorn are internally very similar to factors in R (see ?factor and
?contrasts).

You just have to define your own contrast matrix. Ignoring ambiguous states
for the moment:

contr - matrix(.01, 7, 4, dimnames=list(c(a,c,g,t, -,n,?),
c(a, c, g, t)))
diag(contr) - 1
contr[5:7, ] - 1
contr

contr will look like this
acgt
a 1.00 0.01 0.01 0.01
c 0.01 1.00 0.01 0.01
g 0.01 0.01 1.00 0.01
t 0.01 0.01 0.01 1.00
- 1.00 1.00 1.00 1.00
n 1.00 1.00 1.00 1.00
? 1.00 1.00 1.00 1.00

if you have already an phyDat object x for your nucleotides a command like
this will just change the contrasts:

x - phyDat(as.character(x), type=USER, contrast=contr))

You can use x as input in pml.

A lot of this is actually covered in one of the vignettes ( type
vignette(phangorn-specials) in R). So those who actually RTFM have a
clear advantage ;)
There is an example with gaps as a fifth state. Joe has done this for his
dnapars but not dnaml in phylip.

Storing data like this has a few advantages and it is similar to RAxML I
belive. This way you one use the same ml code for nucleotides, AA, codons,
covarion models or whatever you can think of. Also the information for tips
are stored as integers (not as four doubles) which saves space and one
needs only to subset the contrast matrix, which can save computations.

Hope this helps you get started

Regards,
Klaus

On Thu, Jan 22, 2015 at 1:24 PM, George Shirreff 
georgeshirr...@googlemail.com wrote:

 Hello,

 I was wondering if there is a way to incorporate sequencing error in
 phylogenetic analysis, if I know there is sequencing error but it's
 randomly distributed.

 I imagine this could be as simple as setting the probability of the
 observed base at each tip to 0.99 rather than 1, with a probability of
 0.01/3 rather than 0 for each of the other tips (for a sequencing error of
 1%). I was hoping that something like pml (phangorn) can be modified to do
 that, but I haven't been able to figure it out myself yet.

 Thanks in advance,
 George

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




-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston

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

Re: [R-sig-phylo] Problem loading phytools on Mac

[Klaus Schliep]

Hi Gabriel,

the dependency for the rgl package comes from my phangorn package. You can
download the development version which does not need rgl installed the
following way:
library(devtools)
install_github(KlausVigo/phangorn)

It seems that apple is trying hard not to be a linux system any more and
breaking perfectly working stuff.

Regards,
Klaus


On Mon, Dec 22, 2014 at 4:33 AM, Gabriel Yedid gyedi...@gmail.com wrote:

 Hello list,

 Has anyone else experienced the following error when trying to load
 phytools on a Mac, even though it seems to have been successfully
 installed?  Note that this problem doesn't seem to be with phytools
 itself, but with rgl:

 library(phytools)
  Loading required package: ape
 Loading required package: maps
 Error : .onLoad failed in loadNamespace() for 'rgl', details:
   call: dyn.load(file, DLLpath = DLLpath, ...)
   error: unable to load shared object

 '/Library/Frameworks/R.framework/Versions/3.1/Resources/library/rgl/libs/rgl.so':

 dlopen(/Library/Frameworks/R.framework/Versions/3.1/Resources/library/rgl/libs/rgl.so,
 6): Library not loaded: /opt/X11/lib/libGLU.1.dylib
   Referenced from:

 /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rgl/libs/rgl.so
   Reason: image not found
 Error: package or namespace load failed for 'phytools'



 I haven't experienced this problem using phytools on a PC.  Anybody
 know what gives here?

 cheers,

 Gabe

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




-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston

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

Re: [R-sig-phylo] phangorn - midpoint function

[Klaus Schliep]

Hi Fabricia,

nobody will be able to answer you if you don't provide the data and/or code
you are using. See for example
http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example/5963610#5963610
If you send me the tree, I can try to figure out what is going wrong.

Klaus


On Fri, Nov 21, 2014 at 3:17 PM, Fabricia Nascimento 
nasciment...@yahoo.com.br wrote:

 Hi!

 I am using the function midpoint of package phangorn to root my trees
 (which I generated by simulations).
 A negative branch length is added to my tree after I use the midpoint
 function. The same does not happen if I use FigTree for example.
 Does anyone know why that happens?

 I cannot have a negative branch length on my trees.

 Thanks!
 Fabricia.

 ___
 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-ph...@r-project.org/




-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston

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

Re: [R-sig-phylo] Identify nodes along a pathway from root to a given tip

[Klaus Schliep]

Hello Manabu,
try Ancestors from phangorn.
Regards
Klaus
Am 04.11.2014 08:12 schrieb Manabu Sakamoto manabu.sakam...@gmail.com:

 Dear list,

 I was wondering if there was an existing function to identify all the nodes
 (e.g., phylo$edge[,1]) along a pathway from the root to any given tip.

 I've found Descendants in phangorn, but I only want the node IDs from the
 root leading to a given taxon.

 Many thanks,
 Manabu

 --
 Manabu Sakamoto, PhD
 manabu.sakam...@gmail.com

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

Re: [R-sig-phylo] Network from gene trees

[Klaus Schliep]

Hello Vojt�ch,

have you looked into consensusNet in phangorn?

Regards,
Klaus

Am 03.11.2014 03:10 schrieb Vojt�ch Zeisek vo...@trapa.cz:

 Hello,
 let's say I have many gene trees (all have same labels) in one multiPhylo
 object. The trees are not fully congruent. One of the reasons can be
 hybridizations among taxa. I wonder if I can make from those trees a
 network
 looking, for example, like those from SplitsTree. I looked at functions in
 ape
 and phangorn, but they seem to work on only one phylogenetic tree or on
 distance matrix. I want to get a species tree, but with (possible)
 reticulations. Can I do this?
 Sincerely,
 Vojt�ch

 --
 Vojt�ch Zeisek
 http://trapa.cz/en/

 Department of Botany, Faculty of Science
 Charles University in Prague
 Ben�tsk� 2, Prague, 12801, CZ
 http://botany.natur.cuni.cz/en/

 Institute of Botany, Academy of Science
 Z�mek 1, Pr�honice, 25243, CZ
 http://www.ibot.cas.cz/en/

 Czech Republic

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

Re: [R-sig-phylo] sitewise parsimony scores with parsimony()

[Klaus Schliep]

Dear Rob

pdat has stored the information which site belongs to a site pattern in the
attribute index.
So this should work for you:
parsimony(tree,pdat,site=site)[attr(pdat, index)]

Klaus


On Thu, Oct 30, 2014 at 9:56 AM, Rob Unckless unckl...@cornell.edu wrote:

 Hello,


 I would like to calculate a Fitch parsimony score for each site in an
 alignment using parsimony().  The problem is that the function only returns
 a score for each unique site pattern (as determined by phyDat).  Is there
 anyway (either with phyDat() or parsimony()) to get scores for each site
 regardless of whether or not it is a unique site pattern?


 For example:


  library(ape)
  library(phangorn)
 
  dat-matrix(c(0,1,1,0,0,1,1,0,0,0,0,1,1,1,0,1),ncol=4)
  rownames(dat)-c(A,B,C,D)
  colnames(dat)-c(X1,X2,X3,X4)
  dat
   X1 X2 X3 X4
 A  0  0  0  1
 B  1  1  0  1
 C  1  1  0  0
 D  0  0  1  1
  pdat-phyDat(dat,type=USER,levels=c(0,1))
  pdat
 4 sequences with 4 character and 3 different site patterns.
 The states are 0 1
  tree-rtree(4, rooted = TRUE, tip.label = c(A,B,C,D), br = runif)
  parsimony(tree,pdat,site=site)
 [1] 2 1 1

 There are four positions (X1-X4) but only three site patterns and
 therefore three parsimony scores.

 Thank you.

 Rob





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




-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston

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

Re: [R-sig-phylo] midpoint rooting? how to get the outgroup?

[Klaus Schliep]

Dear Romain,

the small function below returns the out group. I define he outgroup as the
clade with less taxa from the root. If both clades have the same number of
taxa in it  one clades is chosen randomly.

library(phangorn)

getOutgroup - function(tree){
  tree = midpoint(tree) # you probably have done this beforehand
  tmp = Descendants(tree)
  kids = Children(tree, getRoot(tree))
  tmp = tmp[kids]
  x = sapply(tmp, length)
  tree$tip.label[tmp[[which.min(x)]]]
}

# example
trees = rmtree(10, 5)

outGroups - lapply(trees, getOutgroup)
res = sapply(outGroups , function(x, y)any(y %in% x), y=t1)
sum(res) # number of outgroups which contain t1
which(res) # which trees have t1 in the outgroup
which(res  (sapply(outGroups, length)==1)) # which trees have only t1 in
the outgroup

trees = lapply(trees, midpoint)
class(trees) = multiPhylo
plot(trees)

Cheers,
Klaus


PS: midpoint is now part of the phangorn package and much faster for larger
trees.


On Thu, Oct 23, 2014 at 5:57 AM, romain robl...@paca.inra.fr wrote:

 Dear Klauss,

 Following the topic below, I was wondering if there is function to get the
 name of the outgroup after having rooted the tree with the midpoint
 function.
 I can plot the tree and just check it by eyes. However I got  1000 trees
 and I want to count how many times on particular species is placed as
 outgroup.

 Thank you for your help.

 Romain Blanc-Mathieu
 PhD INRA-CNRS-UNS Agrobiotech sophia Antipolis
 06160 Antibes
 FRANCE


 Dear Robin,

 here is a function, that does midpoint rooting, you have to attach the
 phangorn package. It is likely to appear in one of the next ape
 versions.

 Cheers,
 Klaus




 midpoint - function(tree){
  dm = cophenetic(tree)
  tree = unroot(tree)
  rn = max(tree$edge)+1
  maxdm = max(dm)
  ind =  which(dm==maxdm,arr=TRUE)[1,]
  tmproot = Ancestors(tree, ind[1], parent)
  tree = phangorn:::reroot(tree, tmproot)
  edge = tree$edge
  el = tree$edge.length
  children = tree$edge[,2]
  left = match(ind[1], children)
  tmp = Ancestors(tree, ind[2], all)
  tmp= c(ind[2], tmp[-length(tmp)])
  right = match(tmp, children)
  if(el[left]= (maxdm/2)){
   edge = rbind(edge, c(rn, ind[1]))
   edge[left,2] = rn
   el[left] = el[left] - (maxdm/2)
   el = c(el, maxdm/2)
  }
  else{
  sel = cumsum(el[right])
  i = which(sel(maxdm/2))[1]
  edge = rbind(edge, c(rn, tmp[i]))
  edge[right[i],2] = rn
  eltmp =  sel[i] - (maxdm/2)
 #el = c(el, sel[i] - (maxdm/2))
  el = c(el, el[right[i]] - eltmp)
  el[right[i]] = eltmp
  }
  tree$edge.length = el
  tree$edge=edge
  tree$Nnode  = tree$Nnode+1
  phangorn:::reorderPruning(phangorn:::reroot(tree, rn))
 }




 On 9/2/10, Velzen, Robin van Robin.vanVelzen at wur.nl  
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo wrote:
 /  Dear Ape authors and mailing list members,
 //
 //  I am new to R. and to Ape and very much impressed by the
 functionality the
 //  platform offers. I have used Ape to construct distance-based trees
 using the
 //  'read.dna', 'dist.dna' and 'nj' functions. Now I wish to root the
 resulting
 //  Neighbor Joining tree using midpoint rooting, but have not been able
 to find
 //  the right function.
 //
 //  Does anyone know if there is a function for midpoint rooting in Ape or
 //  similar R. package, or if there is a smart way to define a midpoint
 root
 //  outgroup that can be used with the 'root' function in Ape?
 //
 //  Any help or suggestion will be much appreciated!
 //
 //  Thanks,
 //
 //  Robin
 //
 //  Robin van Velzen
 //  PhD student
 //  Biosystematics Group
 //  Wageningen University
 //
 //  Wageningen Campus, Radix building 107, Room W4.Aa.095
 //  Droevendaalsesteeg 1, 6708 PB Wageningen, The Netherlands
 //  PO Box 647, 6700 AP Wageningen, The Netherlands
 //  Tel. +31 (0)317 483425
 //  http://www.bis.wur.nl
 //
 //
 //
 // [[alternative HTML version deleted]]
 //
 //  ___
 //  R-sig-phylo mailing list
 //  R-sig-phylo at r-project.org  
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
 //  https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
 //
 /

 --
 Klaus Schliep
 Universit� Paris 6 (Pierre et Marie Curie)
 9, Quai Saint-Bernard, 75005 Paris


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




-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston

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

Re: [R-sig-phylo] Plotting Posterior Probabilities from MrBayes Trees (An Update)

[Klaus Schliep]

, while leaving the tip.label order
  and the edge lengths the same. Also, larger trees seem to produce
  larger inconsistencies in the tree produced.
 
  So... now for the list, are there any other methods on CRAN for
  obtaining the posterior probabilities from a .con.tre file? Otherwise,
  next I guess I'll be trying phyloch.
 
  Cheers,
  -Dave
 
  --
  David W. Bapst, PhD
  Adjunct Asst. Professor, Geology and Geol. Eng.
  South Dakota School of Mines and Technology
  501 E. St. Joseph
  Rapid City, SD 57701
 
  http://webpages.sdsmt.edu/~dbapst/
  http://cran.r-project.org/web/packages/paleotree/index.html
 
  ___
  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/
 
 
 
 
 
  --
  David W. Bapst, PhD
  Adjunct Asst. Professor, Geology and Geol. Eng.
  South Dakota School of Mines and Technology
  501 E. St. Joseph
  Rapid City, SD 57701
 
  http://webpages.sdsmt.edu/~dbapst/
  http://cran.r-project.org/web/packages/paleotree/index.html



 --
 David W. Bapst, PhD
 Adjunct Asst. Professor, Geology and Geol. Eng.
 South Dakota School of Mines and Technology
 501 E. St. Joseph
 Rapid City, SD 57701

 http://webpages.sdsmt.edu/~dbapst/
 http://cran.r-project.org/web/packages/paleotree/index.html

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




-- 
Klaus Schliep
Postdoctoral Fellow
Revell Lab, University of Massachusetts Boston

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

Re: [R-sig-phylo] Ape package.

[Klaus Schliep]

The tree file is not proper formatted, it is missing a ; in the end.

Regards,
Klaus


On Sun, Apr 20, 2014 at 10:17 PM, Alexey Fomin lesh...@rambler.ru wrote:

  Hi!

 Whether someone work with ape package?
 I have long tree (see attached, it is in newick format, it is ultra
 metric, it is from http://www.onezoom.org/user/tetrapods_tree.js ).
 I need dates of nodes.
 Someone advised me to do:

 1) library(ape)
 2) tree-read.tree(?yourtree.nex?)
 3) br.times-branching.times(tree)

 But I can not not do the second step.

 I tried (after library(ape))
 tree-read.tree(F:\tetrapods_tree.txt)
 But it is gives:
 Error: unexpected input in tree-read.tree(F:\

 What can I do?

 Alexey.

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




-- 
Klaus Schliep

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

Re: [R-sig-phylo] chronos and gammaStat

[Klaus Schliep]

Hello Frabricia,

I assume your RAxML tree is not rooted.
is.rooted(chr)
and use the function root if chr is not rooted.

Regards,
Klaus


On Tue, Apr 8, 2014 at 2:03 PM, Fabricia Nascimento 
nasciment...@yahoo.com.br wrote:

 Hi!

 I am using the gammaStat in ape to analyze some trees I generated using
 simulations and RaxML.


 After using the R code below to convert it to ultrametric tree:
 cal - makeChronosCalib(Tree, age.min=1)
 chr - chronos(Tree, lambda = 0.1, calibration = cal)
 gammaStat(chr)

 I get the following result:
  gammaStat(chr)
 [1] 2.639323
 Warning message:
 In (2:N) * g :
   longer object length is not a multiple of shorter object length

 Can I ignore this warning?

 Thanks very much!
 Fabricia.

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




-- 
Klaus Schliep

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

Re: [R-sig-phylo] chronos and gammaStat

[Klaus Schliep]

You could try
multi2di()
to see if the warning disappears.
It is probably better to do, but I haven't looked into the implementation.


On Tue, Apr 8, 2014 at 2:23 PM, Fabricia Nascimento 
nasciment...@yahoo.com.br wrote:

 Hi!

 Trees are not binary. Do you think this is the problem?
 Should I convert it?

 Thanks!

   --
  *De:* Klaus Schliep klaus.schl...@gmail.com
 *Para:* Fabricia Nascimento nasciment...@yahoo.com.br
 *Cc:* r-sig-phylo@r-project.org r-sig-phylo@r-project.org
 *Enviadas:* Terça-feira, 8 de Abril de 2014 13:20

 *Assunto:* Re: [R-sig-phylo] chronos and gammaStat

 I am guessing. Are the trees all binary, too?

 is.binary.tree(chr)


 On Tue, Apr 8, 2014 at 2:16 PM, Fabricia Nascimento 
 nasciment...@yahoo.com.br wrote:

 Dear Klaus,

 all trees are rooted.
 So I'm not sure why I am getting this warning.

  Thanks!
 Fabricia.

   --
  *De:* Klaus Schliep klaus.schl...@gmail.com
 *Para:* Fabricia Nascimento nasciment...@yahoo.com.br
 *Cc:* r-sig-phylo@r-project.org r-sig-phylo@r-project.org
 *Enviadas:* Terça-feira, 8 de Abril de 2014 13:13
 *Assunto:* Re: [R-sig-phylo] chronos and gammaStat

 Hello Frabricia,

 I assume your RAxML tree is not rooted.
 is.rooted(chr)
 and use the function root if chr is not rooted.

 Regards,
 Klaus


 On Tue, Apr 8, 2014 at 2:03 PM, Fabricia Nascimento 
 nasciment...@yahoo.com.br wrote:

  Hi!
 
  I am using the gammaStat in ape to analyze some trees I generated using
  simulations and RaxML.
 
 
  After using the R code below to convert it to ultrametric tree:
  cal - makeChronosCalib(Tree, age.min=1)
  chr - chronos(Tree, lambda = 0.1, calibration = cal)
  gammaStat(chr)
 
  I get the following result:
   gammaStat(chr)
  [1] 2.639323
  Warning message:
  In (2:N) * g :
   longer object length is not a multiple of shorter object length
 
  Can I ignore this warning?
 
  Thanks very much!
  Fabricia.
 
 [[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/
 



 --
 Klaus Schliep


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





 --
 Klaus Schliep







-- 
Klaus Schliep

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

Re: [R-sig-phylo] Fatal Crash in ancestral.pars in phangorn

[Klaus Schliep]

Thanks, David! I will have a look into it.
Klaus


On Tue, Apr 8, 2014 at 8:43 PM, David Bapst dwba...@gmail.com wrote:

 Hi all,

 I seem to have found a combination of topology and characters that causes
 ancestral.pars in phangorn to hard crash, at least in Windows7 with R v3.03
 and phangorn v1.99-7.

 This issue seems to have been introduced recently; the code works fine in
 phangorn v1.99-5 with R v3.02.

 I have a reproducible example below; as I said, it does hard crash the R
 console, so BEWARE when trying it out. I'm not entirely certain what about
 this particular dataset is causing the crash: the below is a big graptolite
 phylogeny (...as I'm sure some of you guessed). The tree has polytomies,
 but multi2di()-ing the tree doesn't stop the hard crash.

 Cheers,
 -Dave



 library(phangorn)

 char-c(5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,
 5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,3,3,3,3,3,
 3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,
 3,4,3,3,3,3,3,3,3,3,3,3,3,3,3,4,3,3,3,3,3,3,3,3,3,3,3,3,
 3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,3,3,3,3,3,3,4,3,
 3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,5,5,3,3,3,3,3,
 3,3,3,3,3,3,3,3,3,2,1,1,1,3,1,1,3,3,3,3,3,3,3,5,3,3,3,0,
 0,0,0,0,0,0,0,0,0,0,0,0,0,3,3,4,3,3,3,0,0,1,3,3,3,3,0,3,
 3,3,0,3,0,0,0,0,3,4,3,0,0,0,2,0,0,3,0,0,0,0,0,0,0)


 tree-read.tree(file=,text=((t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15,t16,t17,t18,t19,t20,t21,t22,t23,t24,t25,t26,t27,t28,t29,t30,t31,t32,t33,t34,t35,t36,t37,t38,t39,t40,t41,t42,t43,t44,t45,t46,t47),t48),t49,t50),t51),t52,t53),t54,t55,t56,t57,t58),((t59,t60),t61)),t62),(t63,t64)),t65),t66,t67),t68,t69,t70),t71,t72),t73),t74,t75),(t76,t77),t78,t79,t80,t81),(t82,t83),t84,t85,t86),t87,t88),t89),t90,t91,t92,t93,t94),t95,t96,t97,t98,t99,t100,t101,t102),t103),t104),t105),(t106,t107),t108,t109,t110,t111),t112),t113),(t114,t115)),t116),t117),t118),t119,t120),t121),((t122,t123),t124)),(t125,t126)),t127),t128,t129),t130),t131)),t132),t133),t134),t135,t136,t137,t138),t139,t140,t141,t142),t143),t144),(t145,t146),t147),t148),t149),t150,t151))),t152,t153),t154),t155),t156),t157),t158,t159),t160),t161,t162,t163,t164,t165),t166,t167,t168)),t169),((t170,t171,t172),t173)),(((t174,t175,t176),t177),(t178,t179)),t!
 
180),(((t181,t182),t183),t184,t185),(t186,t187,t188,t189),t190,t191),t192,t193,t194,t195,t196,t197,t198,t199,t200,t201,t202,t203,t204,t205),(((t206,t207,t208,t209,t210,t211),t212,t213,t214),(t215,t216,t217,t218,t219,t220,t221),(t222,t223,t224),t225,t226,t227,t228,t229,t230,t231,t232,t233,t234)),t235,t236,t237,t238,t239,t240,t241,t242,t243,t244,t245);)

 char-matrix(char,,1)
 rownames(char)-tree$tip.label
 char-phyDat(char,type=USER,levels=sort(unique(char)))

 anc-ancestral.pars(tree,char,type=MPR)  #here's the crash!



 --
 David W. Bapst, PhD
 Adjunct Asst. Professor, Geology and Geol. Eng.
 South Dakota School of Mines and Technology
 501 E. St. Joseph
 Rapid City, SD 57701

 http://webpages.sdsmt.edu/~dbapst/
 http://cran.r-project.org/web/packages/paleotree/index.html




-- 
Klaus Schliep

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

Re: [R-sig-phylo] Generate all possible trees for n species

[Klaus Schliep]

There is a function allTrees in phangorn, which computes all binary trees.

Cheers,
Klaus


On Sat, Mar 8, 2014 at 1:40 PM, Augusto Ribas ribas@gmail.com wrote:

 Hello.

 I'm trying to do a function to generate all possible trees for N species.
 I think i almost got it, but there is something that i'm missing.

 ##
 todas_arvores-function(entrada) {
 #Transform the vector of names in a list with one element
 saida-list(entrada)
 #In the fisrt iteration, the first element of this list will have
 length==entrada,
 #we will iterate until it is one full string, just one element
 while(length(saida[[1]])!=1) {
 #auxiliar list to receive the sub-trees
 aux-list()
 #the first iteration this is one, but increase within each
 iteration
 for(i in 1:length(saida)) {
 #another temporarily container for sub-trees, here to be able
 to iterate though all the saida length
 temporario-list()
 #Do all possible 2 combinations, unite 2 species, or one
 species and a sub-tree(now a greater string)
 #(a,b) for example, after the first iteration and reserve here
 combinacoes-combn(saida[[i]],2)
 #now for all combinations i glue the combination with paste,
 like (a,b), and generate a vector with what is
 #still missing, look the %in% check
 for(j in 1:ncol(combinacoes)) {


 temporario[[j]]-c(saida[[i]][!saida[[i]]%in%combinacoes[,j]],paste((,combinacoes[1,j],,,combinacoes[2,j],),sep=))
 }
 #for each element os saida, i store the results here
 aux-c(aux,temporario)
 }
 #at the end, i store everything at saida again, now we are ready to
 test the while condition again, but now i
 #hope that there will be a vector of n-1 elements from the when we
 started this iteration.
 saida-aux
 }
 #just to produce a better output, i unlist and add the ; in the end, to
 use read.tree from ape
 saida-unlist(saida)
 saida-sapply(saida,function(arvore)
 {paste(arvore,;,sep=)},USE.NAMES = FALSE)
 return(saida)
 }


 #example of use
 arvores-todas_arvores(paste(sp,1:4))
 arvores

 library(ape)

 #example of wrong output
 dist.topo(read.tree(text=arvores[1]),read.tree(text=arvores[16]))


 #just a plot to see what i'm trying
 raiz-sqrt(length(arvores))
 par(mar=c(1,1,1,2))

 layout(matrix(1:c(floor(raiz)*ceiling(raiz)),nrow=floor(raiz),ncol=ceiling(raiz),byrow=T))

 for(i in 1:length(arvores)) {
 plot(read.tree(text=arvores[i]),cex=1)
 title(paste(Tree,i))
 }

 ##


 as i see, here

 for(j in 1:ncol(combinacoes)) {


 temporario[[j]]-c(saida[[i]][!saida[[i]]%in%combinacoes[,j]],paste((,combinacoes[1,j],,,combinacoes[2,j],),sep=))
 }

 i would need to test if the sub-tree(elements of the vector) is already
 present in the aux list, to evade reintroduce copy of the same tree, but
 i don't know how.
 something like, there will be a moment, for 4 species that there will be a
 element in aux like
 (sp 1,sp 2)  (sp 3,sp 4)

 this time i dont want to put this one on temporario
 (sp 3,sp 4) (sp 1,sp 2)

 because there is already an element like this on aux, just the inverse
 order, as i keep both i produce 2 equal trees in the end.


 I don't know if maybe there is a simple algorithm to do this, something
 like apply a modification to a starting tree, any suggestion?

 Another question is, i need to produce the species catalan number of trees
 right? Is this how to predict the number of possible trees?

 Best Wishes, Augusto Ribas

 --
 Grato
 Augusto C. A. Ribas

 Site Pessoal: http://recologia.com.br/ http://augustoribas.heliohost.org
 Github: https://github.com/Squiercg
 Lattes: http://lattes.cnpq.br/7355685961127056

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




-- 
Klaus Schliep
Phylogenomics Lab at the University of Vigo, Spain
http://darwin.uvigo.es/kschliep/

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

Re: [R-sig-phylo] Midpoint rooting routine?

[Klaus Schliep]

Hi all,
I just uploaded a new version of phangorn (1.99-7), which does not depend
on the rgl package. This should allow phangorn to work on all platforms
again. It seems midpoint (phangorn) and midpoint.root (phytools) give the
same result, which is reassuring. midpoint seems to be faster for larger
trees, let say more than 1000 tips.
I also started a small side project (
http://klash.shinyapps.io/shinyTreeViewer/) which is an interactive
web-based tree viewer using mainly the ape, phangorn  shiny package. There
is a web-server version from
http://klash.shinyapps.io/shinyTreeViewer/hosted kindly by RStudio.
The github version is more up to date. The tree
viewer generates also some R-code to either use it as starting point to
generate more complex graphics and/or for reproducible trees. Feedback is
welcome.
Regards,
Klaus



On Fri, Jan 31, 2014 at 1:11 AM, Todd Oakley
todd.oak...@lifesci.ucsb.eduwrote:

 All,
 As an update: the script that Liam wrote for midpoint rooting seems to
 have worked well for us. I did not numerically test the rooting results,
 but they pass the eye test (ie they look right to me). I serially rooted
 about 100 trees with 50-250 OTU's each, and it was fast, so it was
 efficient in our hands.

 Liam has described his process and function on his blog, here:

 http://blog.phytools.org/2014/01/function-for-midpoint-
 rooting.html?spref=fb



 My application for this is a tool called tab2trees, implemented in Osiris
 for Galaxy. It takes a list of newick tree descriptions and visualizes all
 the trees in a pdf book, with one tree per page. Since I often use gene
 trees, midpoint rooting is a sensible way to present the trees. Here is a
 description of tab2trees (pre-midpoint rooting):

 http://osiris-phylogenetics.blogspot.com/2012/09/tab2trees.html


 Best wishes,
 Todd Oakley, UCSB






 On 1/30/2014 3:57 PM, Joe Felsenstein wrote:

 Liam said:

  Joe's correct that this is not yet in our project to create an R
 interface
 for PHYLIP (Rphylip, https://github.com/liamrevell/Rphylip). As soon as
 I
 get a chance to work on this project today, I will add it.

 That may be hard to do, as Retree is interactive and
 has two levels of menus.

 Joe
 
 Joe Felsenstein j...@gs.washington.edu
   Department of Genome Sciences and Department of Biology,
   University of Washington, Box 355065, Seattle, WA 98195-5065 USA


 ___
 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-ph...@r-project.org/




-- 
Klaus Schliep
Phylogenomics Lab at the University of Vigo, Spain
http://darwin.uvigo.es/kschliep/

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

Re: [R-sig-phylo] Run time for pglmm.fit in picante

[Klaus Schliep]

Dear Dan,

it is best to do some profiling of the code, there are useful functions
Rprof() and summaryRprof(). In most cases your covariance matrix and
dependent variables will be highly structured and in most cases sparse. You
can use the Matrix package to handle sparse matrices. This will save a lot
of computations and can speed up computations several orders of magnitude
if you are lucky.
Along the same lines the function eigen checks if the matrix is symmetric
or not. This should be known beforehand or can be done out side the
optimize function.

Cheers,
Klaus


On Fri, Sep 6, 2013 at 6:44 PM, Dan Larkin dlar...@chicagobotanic.orgwrote:

 Hello,

 I'm encountering long running times using the pglmm.fit function in
 Picante. Ives  Helmus warn that the estimation procedure can be slow for
 large datasets so this is not surprising. But my question is whether the
 times I am experiencing sound like normal slow or weird slow? And does
 anyone have suggestions for tweaks that might speed things up?

 Information:
 - System: Windows 7, 16 GB RAM, R 3.0.1 (64-bit)
 - Running IH's Model I: Is there a phylogenetic pattern in community
 structure?
 - Dataset: 34 sites, 236 angiosperm species

 Run times:
 - Using the pglmm sample code (simulated data for 30 sites w/ 16 spp from
 a balanced tree, exitcountermax = 30 [iterations for fixed and random
 effects estimation]):
 # user  system elapsed
 # 12.341.56   14.18

 - My data arbitrarily trimmed down to same (30 sites, 16 spp, same
 exitcountermax)
 # user  system elapsed
 # 32.873.13   36.05

 - My data trimmed to 30 sites, 32 spp, same exitcountermax
 # user  system elapsed
 # 267.45   15.19  282.90

 - When I try to fit the model using the full dataset (34 sites, 236 spp),
 each iteration (1 exitcountermax) takes ~4-6 hrs.

 Do these times more-or-less fit expectations? As far as I can tell, it's
 where pglmm.fit uses the 'optim' and 'optimize' functions to estimate
 variances of the random effects w/ REML that is slow.

 Ultimately, I'd like to use some of the more complex models developed by
 IH (2011) and will likely need high iterations to achieve convergence.
 Does anyone have ideas for ways I could modify the pglmm.fit function to
 speed the estimation process?

 Thanks in advance for any suggestions. And thanks to Tony and Matt for
 developing these models and making them available in R.

 Dan

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




-- 
Klaus Schliep
Phylogenomics Lab at the University of Vigo, Spain
http://darwin.uvigo.es/kschliep/

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

Re: [R-sig-phylo] compar.gee crashing R

[Klaus Schliep]

Dear Sereina,

you should give an example to reproduce the error and when possible the
output to give others the opportunity to locate your problem. Also run
traceback(), if just the function crashes.

Regards,
Klaus


On Sat, Aug 24, 2013 at 3:08 PM, Sereina Graber sereina.gra...@gmx.chwrote:

 Dear r-sig-phylo users,

 I am running a phylogenetic logistic regression using the function
 compar.gee() in the  package ape. Sometimes, when running such a
 phylogenetic logistic regression, RStudio crashes. However, it never
 happened to me when using compar.gee for modelling a normal continuous
 trait. It seems it has nothing to do with the strength of the correlation
 coeffiicent, so nothing like a perfect fit problem. Did anyone of you have
 the same problem with crashing of R? what might be the problem here?
 converging problems?

 Thank you  cheers,
 Sereina

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




-- 
Klaus Schliep
Phylogenomics Lab at the University of Vigo, Spain
http://darwin.uvigo.es/kschliep/

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

Re: [R-sig-phylo] Collapsing branches with low bootstrap values

[Klaus Schliep]

On Thu, Aug 22, 2013 at 5:09 PM, Naxerova, Kamila
naxer...@fas.harvard.eduwrote:

 Dear list,

 is there a function that will take a tree and the output of boot.phylo()
 and collapse branches below a chosen bootstrap value cutoff?

 Many thanks.
 Kamila

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




-- 
Klaus Schliep
Phylogenomics Lab at the University of Vigo, Spain
http://darwin.uvigo.es/kschliep/

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

[R-sig-phylo] Fwd: Collapsing branches with low bootstrap values

[Klaus Schliep]

-- Forwarded message --
From: Klaus Schliep klaus.schl...@gmail.com
Date: Thu, Aug 22, 2013 at 6:12 PM
Subject: Re: [R-sig-phylo] Collapsing branches with low bootstrap values
To: Naxerova, Kamila naxer...@fas.harvard.edu


Hi Kamila,

try function pruneTree in phangorn.
data(woodmouse)
f - function(x) nj(dist.dna(x))
tr - f(woodmouse)
X = boot.phylo(tr, woodmouse, f, trees = TRUE)
tree = plotBS(tr, X$trees)
tree2 = pruneTree(tree, 75)
par(mfrow=c(2,1))
plot(tree, show.node=TRUE)
plot(tree2, show.node=TRUE)

Cheers,
Klaus


On Thu, Aug 22, 2013 at 5:09 PM, Naxerova, Kamila
naxer...@fas.harvard.eduwrote:

 Dear list,

 is there a function that will take a tree and the output of boot.phylo()
 and collapse branches below a chosen bootstrap value cutoff?

 Many thanks.
 Kamila

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




-- 
Klaus Schliep
Phylogenomics Lab at the University of Vigo, Spain
http://darwin.uvigo.es/kschliep/




-- 
Klaus Schliep
Phylogenomics Lab at the University of Vigo, Spain
http://darwin.uvigo.es/kschliep/

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

Re: [R-sig-phylo] getting a list of clusters with their associated bootstrap support

[Klaus Schliep]

Hello,


On Tue, Jun 11, 2013 at 5:53 PM, Douda Bensasson 
douda.bensas...@manchester.ac.uk wrote:


 I have generated neighbour joining distance trees in ape, and can plot or
 write them with their associated bootstrap support. However for various
 downstream analyses I would also like to get a list of clusters with their
 associated bootstrap support, in a text form that is more easy to parse
 than a newick format tree.

 I thought that prop.part would do this, but the groupings and frequencies
 I get using prop.part differ from what I see when I plot trees in R (or
 write them in Newick format).

 This is what I have tried so far:

 library(ape)
 data(woodmouse)
 f - function(x) nj(dist.dna(x))
 tr - f(woodmouse)
 set.seed(1)
 bs-boot.phylo(tr, woodmouse, f, quiet = TRUE, trees=TRUE)
 plot(tr)
 nodelabels(bs$BP)
 tr$node.label-bs$BP
 write.tree(tr,njwithbsWM.tree)
 prop.part(bs$trees,check.labels=TRUE)

 Try:
pp - prop.part(bs$trees,check.labels=TRUE) # contains a list of vectors
with the tips descending from each node
pp2 - postprocess.prop.part(pp) # contains bipartitions
prop.clades(tr, part = pp) # bs values for the tree tr

Regards,
Klaus





 .. and the list of descendent tips does not correspond to the frequencies
 seen in the node labels of the trees I plot.

 Any help would be much appreciated!

 Douda
 ___
 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/




-- 
Klaus Schliep
Phylogenomics Lab at the University of Vigo, Spain
http://darwin.uvigo.es/kschliep/

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

Re: [R-sig-phylo] Incongruence

[Klaus Schliep]

Hi Marcelo,

unfortunately incongruence is not well defined. Here are several
possibilities to access incongruence,
with phangorn/ ape. I would first check if the base frequencies are similar
in each partitions.
There are a few exploratory tools for displaying incongruence like split
networks or the densiTree
to compare trees from different partitions.
More formally you can use a SH-test (SH.test) to test if trees or models
estimated for partitions differ.

Regards,
Klaus



On Tue, May 14, 2013 at 6:12 AM, Marcelo Reginato
reginato...@yahoo.com.brwrote:

 Hi,
 I'd like to know if anyone is aware of any tool in R to assess
 incongruence among dna partitions.
 Best,
 Marcelo

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




-- 
Klaus Schliep
Phylogenomics Lab at the University of Vigo, Spain
http://darwin.uvigo.es/kschliep/

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

Re: [R-sig-phylo] Very simple question on fitted values and scatterplots

[Klaus Schliep]

Hi Xavier

On 5/2/13, Xavier Prudent prudentxav...@gmail.com wrote:
 Dear all,

 I went through this tutorial on the fit of a linear model using Independent
 Contrasts
 http://nunn.rc.fas.harvard.edu/groups/pica/wiki/c57f6/

 What is the meaning of fitted values in the scatter plot (bottom,
 x-axis)?

 Fitting a linear model with a constraint to the origin, so y = a * x, the
 only fitted parameter is a.

The linear model is
 y = a * x + e, with e is N(0, sigma^2)
The fitted values are the yhat:
yhat = a * x

so it makes sense to have several values.

Regards,
Klaus

 Which does not make sense as we would not have one value per node, and we
 would not expect these points to be spread for a good fit.

 Thanks in advance,

 regards,
 Xavier


 --
 *---
 Xavier Prudent
 *
 *
 Computational biology and evolutionary genomics
 *
 *
 *
 *Guest scientist at the Max-Planck-Institut für Physik komplexer Systeme*
 *(MPI-PKS)*
 *Noethnitzer Str. 38*
 *01187 Dresden
 *
 *
 *
 *Max Planck-Institute for Molecular Cell Biology and Genetics*
 *
 (MPI-CBG)
 *
 *
 Pfotenhauerstraße 108
 *
 *
 01307 Dresden
 *
 *

 *
 *
 Phone: +49 351 210-2621
 *
 *Mail: prudent [ at ] mpi-cbg.de
 **---*

   [[alternative HTML version deleted]]




-- 
Klaus Schliep
Phylogenomics Lab at the University of Vigo, Spain

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

Re: [R-sig-phylo] read.nexus.data parser

[Klaus Schliep]

Hi all,

I had a short look in the code and found some bits to speed the
read.nexus.data function up. I added Emmanuel on the list so he may
can put it into the next ape release if it does work.
Generally I agree with Johan that if speed matters fasta files are the
way to go. Nexus files are ugly to parse and contain many
inconsistencies, like parameters inside comments [].

Regards,
Klaus


On 4/9/13, Johan Nylander johan.nylan...@abc.se wrote:
 Dear All,

 Just to avoid confusion, the readNexus function is in the phylobase
 package. And as Ben pointed out, other packages have their own functions
 for reading the data part from a nexus-formatted file, see e.g., read.nex
 in phyloch.

 On a related note, I wrote read.nexus.data as a temporary, crude parsing
 function while waiting for the phylobase project to take off (phylobase
 uses NCL by Lewis  Holder - _the_ nexus parser), so expect
 read.nexus.data to have it's limitations.

 Furthermore, if speed is the concern, it would perhaps be preferable to
 first convert the Nexus data to Fasta, and then use one of the many
 fast(er) parsers implemented in numerous R packages.

 Cheers
 Johan


 On 04/07/2013 02:59 PM, Ben Bolker wrote: On 13-04-05 01:29 PM, Jessica
 Sabo wrote:
 Hi All,

 I am wondering if there is anyway to increase the speed of the
 read.nexus.data parser. Or if there is an alternative that is a
 faster nexus file data parser.

 THanks, Jess


I don't know if it's faster or not, but there is ?readNexus in the
 'ape' package.  Also see library(sos); findFn(read {nexus format})

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


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



-- 
Klaus Schliep
Phylogenomics Lab at the University of Vigo, Spain

read.nexus.data - function (file) 
{
  find.ntax - function(x) {
for (i in 1:NROW(x)) {
  if (any(f - grep(\\bntax, x[i], ignore.case = TRUE))) {
ntax - as.numeric(sub((.+?)(ntax\\s*\\=\\s*)(\\d+)(.+), 
   \\3, x[i], perl = TRUE, ignore.case = TRUE))
break
  }
}
ntax
  }
  find.nchar - function(x) {
for (i in 1:NROW(x)) {
  if (any(f - grep(\\bnchar, x[i], ignore.case = TRUE))) {
nchar - as.numeric(sub((.+?)(nchar\\s*\\=\\s*)(\\d+)(.+), 
\\3, x[i], perl = TRUE, ignore.case = TRUE))
break
  }
}
nchar
  }
  find.matrix.line - function(x) {
for (i in 1:NROW(x)) {
  if (any(f - grep(\\bmatrix\\b, x[i], ignore.case = TRUE))) {
matrix.line - as.numeric(i)
break
  }
}
matrix.line
  }
  trim.whitespace - function(x) {
gsub(\\s+, , x)
  }
  trim.semicolon - function(x) {
gsub(;, , x)
  }
  X - scan(file = file, what = character(), sep = \n, quiet = TRUE, 
comment.char = [, strip.white = TRUE)
  ntax - find.ntax(X)
  nchar - find.nchar(X)
  matrix.line - find.matrix.line(X)
  start.reading - matrix.line + 1
  Obj - vector(list, ntax)
  for(i in 1:ntax)Obj[[i]] = rep(NA, nchar)
  
  i - 1
  pos - 0
  tot.nchar - 0
  tot.ntax - 0
  for (j in start.reading:NROW(X)) {
Xj - trim.semicolon(X[j])
if (Xj == ) {
  break
}
if (any(jtmp - grep(\\bend\\b, X[j], perl = TRUE, ignore.case = TRUE))) {
  break
}
ts - unlist(strsplit(Xj, (?=\\S)(\\s+)(?=\\S), perl = TRUE))
#browser()
if (length(ts)  2) {
  stop(nexus parser does not handle spaces in sequences or taxon names (ts2))
}
if (length(ts) != 2) {
  stop(nexus parser failed to read the sequences (ts!=2))
}
Seq - trim.whitespace(ts[2])
Name - trim.whitespace(ts[1])
nAME - paste(c(\\b, Name, \\b), collapse = )

if (any(l - grep(nAME, names(Obj {
  tsp - strsplit(Seq, NULL)[[1]]
  
  Obj[[l]][pos + c(1:length(tsp))] - tsp
  chars.done - length(tsp)   
  
}
else {
  names(Obj)[i] - Name
  tsp - strsplit(Seq, NULL)[[1]]
  
  Obj[[i]][pos + c(1:length(tsp))] - tsp
  chars.done - length(tsp)  
  
}
tot.ntax - tot.ntax + 1
if (tot.ntax == ntax) {
  i - 1
  tot.ntax - 0
  tot.nchar - tot.nchar + chars.done
  if (tot.nchar == nchar * ntax) {
print(ntot was more than nchar*ntax)
break
  }
  pos - tot.nchar
}
else {
  i - i + 1
}
  }
  if (tot.ntax != 0) {
cat(ntax:, ntax, differ from actual number of taxa in file?\n)
stop(nexus parser did not read names correctly (tot.ntax!=0))
  }
  for (i in 1:length(Obj)) {
if (length(Obj[[i]]) != nchar) {
  cat(names(Obj[i]), has, length(Obj[[i]]), characters\n

Re: [R-sig-phylo] problem with drop.tip

[Klaus Schliep]

Dear John,

can you please be a bit more specific with your error message. It is
always good to have a reproducible example, e.g. adding a tree which
where drop.tip fails, and to run traceback() just after the error to
get more information where the error occurred. It is also useful to
add information of the version of ape and your operating system.

Regards,
Klaus




On 2/28/13, john d dobzhan...@gmail.com wrote:
 Dear all,

 I'm trying to prune a set of 1000 post-burnin trees to include only a
 subset of taxa. Unfortunately the tree is too big to send to the list,
 but if it is really necessary I'll figure out a way to do it.

 tr is my tree and taxa is my list of selected terminals.

 for(i in 1:1000){
   write.tree(drop.tip(tr[[i]],tr[[i]]$tip.label[-match(taxa,
 tr[[i]]$tip.label)]), file=result.tre, append=TRUE)
 }

 If I run that code, it works for some trees, but not for others, for
 which I got the message Error in kids[[parent[i]]] : subscript out of
 bounds.

 Any suggestions?

 John

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



-- 
Klaus Schliep
Phylogenomics Lab at the University of Vigo, Spain

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

Re: [R-sig-phylo] problem with drop.tip

[Klaus Schliep]

Hi John,
It seems a problem occurs within write.tree and not with the pruning.
So try prune the trees first and then write them out.

tr2=vector(list, 1000)
for(i in 1:1000){
  tr2[[i]] - drop.tip(tr[[i]],tr[[i]]$tip.label[-match(taxa,
tr[[i]]$tip.label)])
}
class(tr2) = multiPhylo
plot(tr2)
for(i in 1:1000){print(i);write.tree(tr2[[i]])} # may helps find you
the trees which fail
write.tree(tr2, file=result.tre)

Cheers,
Klaus


On 2/28/13, Klaus Schliep klaus.schl...@gmail.com wrote:
 Dear John,

 can you please be a bit more specific with your error message. It is
 always good to have a reproducible example, e.g. adding a tree which
 where drop.tip fails, and to run traceback() just after the error to
 get more information where the error occurred. It is also useful to
 add information of the version of ape and your operating system.

 Regards,
 Klaus




 On 2/28/13, john d dobzhan...@gmail.com wrote:
 Dear all,

 I'm trying to prune a set of 1000 post-burnin trees to include only a
 subset of taxa. Unfortunately the tree is too big to send to the list,
 but if it is really necessary I'll figure out a way to do it.

 tr is my tree and taxa is my list of selected terminals.

 for(i in 1:1000){
   write.tree(drop.tip(tr[[i]],tr[[i]]$tip.label[-match(taxa,
 tr[[i]]$tip.label)]), file=result.tre, append=TRUE)
 }

 If I run that code, it works for some trees, but not for others, for
 which I got the message Error in kids[[parent[i]]] : subscript out of
 bounds.

 Any suggestions?

 John

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



 --
 Klaus Schliep
 Phylogenomics Lab at the University of Vigo, Spain



-- 
Klaus Schliep
Phylogenomics Lab at the University of Vigo, Spain

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

Re: [R-sig-phylo] Does anyone know about a function to remove redundancy from a database?

[Klaus Schliep]

Hi Luiz,

phangorn has a generic unique function for phyDat objects.
There may be more elegant ways, but this code should work:
library(phangorn)
y = as.DNAbin(unique(as.phyDat(x)))

Cheers,
Klaus


On 2/19/13, Luiz Max Carvalho luizepidemiolo...@gmail.com wrote:
 I've been looking for a function that could scan a DNAbin object and return
 just the unique sequences. Does anyone know about a function for that?

 Thanks in advance,

 --
 Luiz Max Fagundes de Carvalho, IC, Programa de Computação Científica
 (PROCC), Fundação Oswaldo Cruz, Rio de Janeiro, Brasil.
 https://www.researchgate.net/profile/Luiz_Fagundes_de_Carvalho/?ev=prf_info

   [[alternative HTML version deleted]]




-- 
Klaus Schliep
Phylogenomics Lab at the University of Vigo, Spain

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

Re: [R-sig-phylo] issue with cophenetic.phylo()

[Klaus Schliep]

Hi Kelly,

where is the problem? This happens if you work with floating point numbers.
There is a nice link on the developer.r-project on the topic
www.validlab.com/goldberg/paper.pdf

The value is on my machine actually the same as
.Machine$double.eps
[1] 2.220446e-16
which is the smallest positive floating-point number ‘x’ such that ‘1 + x != 1’.

Regards,
Klaus


On 12/7/12, Kelly Burkett kelly.burk...@mail.mcgill.ca wrote:
 Hi Everyone,

 I have been using the cophenetic.phylo() function and I have been
 having some issues with pairs of tips having distances that are
 different than what I would expect.

 The following code illustrates the issue (I have set the seed in
 order to reproduce these results). I generate a tree with
 10 tips. There should be 9 unique non-zero distances in the matrix
 returned by cophenetic.phylo(). Instead there are 10 because
 1.22313194 is duplicated:

 set.seed(1723)
 tree=rcoal(10)
 x=cophenetic.phylo(tree)
 y=sort(unique(as.vector(x)))

 y:
 [1] 0. 0.04371057 0.16334724 0.22292128 0.31957449 0.33656861
 [7] 0.35227827 0.84410552 1.22313194 1.22313194 2.14497976

 As an example, t9 should have the same distance to t8 and t2.
 x[t9,]:
t1t3t6t9t7t2t8
 t5
 0.3522783 0.3365686 0.3365686 0.000 1.2231319 1.2231319 1.2231319
 2.1449798
   t10t4
 2.1449798 2.1449798

 But they are different:
 x[t9,t8]==x[t9,t2]
 [1] FALSE

 They are off by a very small amount:
 x[t9,t8]-x[t9,t2]
 [1] -2.220446e-16

 I am posting this in case it is of interest to others. As a quick fix,
 I am now rounding x.

 Thanks,
 Kelly

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



-- 
Klaus Schliep
Phylogenomics Lab at the University of Vigo, Spain

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

Re: [R-sig-phylo] estimating mutation rates

[Klaus Schliep]

Hi John,

optim.pml has an parameter optRooted and you have to set it to TRUE.
If you give the function an ultrametric tree it will optimise the edge
lengths. Unfortunately it is not possible to use tree rearrangements
with this option so far.

Regards,
Klaus

On 11/9/12, john d dobzhan...@gmail.com wrote:
 Thanks, guys. I just have a follow-up question: Does anyone know of a
 way to optimize BLs using optim.pml() while enforcing a molecular
 clock?

 Thanks!

 John

 On Thu, Nov 8, 2012 at 6:06 PM, Rob Lanfear rob.lanf...@gmail.com wrote:
 Hi John,

 Brian's suggestion is a good one.

 Another option would be to use any of a number of molecular dating
 programs
 (not strictly an R response, sorry) like R8S (for ML) or BEAST (for
 Bayesian). You can then specify your node dates and topology as
 constraints,
 and use the molecular data appropriately. The advantage of these
 approaches
 is that you'll have full control over the models of rate change, and over
 the treatment of uncertainty in the dates and rates (depending on which
 method you choose).

 However, if you really do believe your dates without question, and if
 your
 molecular dataset has sufficient power for rate estimates (i.e. at least
 a
 handful of substitutions estimated on each branch), I would just do as
 Brian
 suggested.

 Rob


 On 9 November 2012 06:35, Brian O'Meara bome...@utk.edu wrote:

 If you want the branch-specific substitution rate, you could estimate
 the
 branch lengths on a topology constrained to match the chronogram, then
 divide each of the molecular tree's branch lengths by the corresponding
 chronogram branch lengths to get these rates. If you want a global rate,
 you might constrain the molecular tree to have the same proportional
 branch
 lengths as the chronogram. Look at pml() in phangorn for calculating the
 likelihoods.

 Best,
 Brian

 ___
 Brian O'Meara
 Assistant Professor
 Dept. of Ecology  Evolutionary Biology
 U. of Tennessee, Knoxville
 http://www.brianomeara.info

 Students wanted: Applications due Dec. 15, annually
 Postdoc collaborators wanted: Check NIMBioS' website
 Calendar: http://www.brianomeara.info/calendars/omeara


 On Thu, Nov 8, 2012 at 12:49 PM, john d dobzhan...@gmail.com wrote:

  Hi there,
 
  I have a well-supported phylogeny (with branch lengths proportional to
  absolute time) and I'm interested in estimating the mutation rate of
  an independently-obtained molecular dataset. Any suggestions?
 
  thanks!
 
  John
 
  ___
  R-sig-phylo mailing list
  R-sig-phylo@r-project.org
  https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
 

 [[alternative HTML version deleted]]


 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo




 --
 Rob Lanfear
 Research Fellow,
 Ecology, Evolution, and Genetics,
 Research School of Biology,
 Australian National University

 www.robertlanfear.com

 July-September
 National Evolutionary Synthesis Center,
 2024 W. Main Street,
 Suite A200,
 Durham, NC 27705-4667,
 USA

 phone: +1 919 668 4592


 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo



-- 
Klaus Schliep
Phylogenomics Lab at the University of Vigo, Spain

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] problem with boot.phylo() and sample order randomisation

[Klaus Schliep]

Hi Alastair,

I discovered two problems with the analysis. One is your data, which
contains too many transitions and so you lost most of the information
about your the process which generated the data.
There is also a problem with prop.clades, which does not take care
enough for different orderings of the tip labels of the input trees.


Here is a small function, you will need:


checkLabels - function(tree, tip){
ind - match(tip, tree$tip.label)
tree$tip.label - tree$tip.label[ind]
ind2 - match(1:length(ind), tree$edge[, 2])
tree$edge[ind2, 2] - order(ind)
tree
}

I edited the text below.

Regards,
Klaus

On 7/8/12, Alastair Potts pott...@gmail.com wrote:
 Hi all,

 I was wondering whether anyone could offer advice on the following problem.

 When samples are closely related (say identical sequences) are
 analysed using the boot.phylo() function, the bootstrap support
 between these samples will be 100% explicitly because the sample order
 in boot.phylo() is not randomised. This means that closely related
 samples are always built with the same branching order even though the
 branch length may be 0. As the branching order is unchanged these
 branches are given high bootstrap values.

Is this really a problem? Maybe you just use di2multi to collapses the
corresponding dichotomies into a multichotomy.


 In previous version of the ape library, including the sample()
 function (see below) was sufficient to overcome this problem. However,
 when running some old analyses I am now getting very different
 results. Including the sample() function now leads to no bootstrap
 support for any branches.

 I have been playing with all combinations of the parameter settings
 that I can find for the boot.phylo() function, but can't seem to find
 a solution.

 I have included sample code below which I believe accurately
 represents the problem I am observing in my real datasets. Any
 suggestions or comments would be greatly appreciated.

 Thank you for your time,
 Alastair Potts

 library(ape)
 library(phangorn)

 windows()
 par(mfrow=c(1,2))

 # Create random tree, set some branch lengths to 0, and simulate DNA
 set.seed(1)
 tr - rtree(20)
 plot(tr)
 #nodelabels()
 tr$edge.length[c(30:33,15:25)] -0
 plot(tr)

 # The current figure shows the random tree before (left) and after (right)
 # shortening of branch lengths

 # Simulate DNA on the tree with shortened branch lengths
 set.seed(1)
 dna - as.DNAbin(simSeq(tr))


# If you try
cophenetic(tr)
# you see that many distances are 4, and so you expect 4 transitions
per site between this
# sequences.
dist.dna(dna, raw)
# contains many distances close to .75. If you sample these it is
likely you generate
# many NaNs later on in
dist.dna(dna, pairwise.deletion = T)
# or having a distance close to 750 in this case
dist.dna(dna,model=N, pairwise.deletion = T)
# For this reason:
tr$edge.length - tr$edge.length / 5
set.seed(1)
dna - as.DNAbin(simSeq(tr))




 ### Perform NJ and bootstrap without randomising the sample order
 # This incorrectly calculates high bootstrap values for short samples
 on short branches
 f - function(x) nj(dist.dna(x,model=N, pairwise.deletion = T))
 tr1 - f(dna)
 bb - boot.phylo(tr1, dna, f ,B=100, trees=TRUE, rooted=T)
 tr1$node.labels - prop.clades(tr1, bb$trees,rooted=F)

 plot(tr1);title(Without sample order randomisation)
 nodelabels(tr1$node.label)

 # Note that samples on incredibly short branch lengths have very high
 # bootstrap values


 ### Perform NJ and bootstrap WITH randomising the sample order
 # This used to work, but now does something weird. Using sample() to
 randomise the
 # input order. See below.
 f2 - function(x) nj(dist.dna(x[sample(1:nrow(x)),],model=N,
 pairwise.deletion = T))
 tr2 - f2(dna)
 bb2 - boot.phylo(tr2, dna, f2 ,B=100, trees=TRUE, rooted=T)
#  tr2$node.labels - prop.clades(tr2, bb2$trees,rooted=T)
# Here comes the fix
btr2 - .compressTipLabel(bb2$trees)
tr2 - checkLabels(tr2, attr(btr2, TipLabel))
btr2 - .uncompressTipLabel(btr2)

tr2$node.labels - prop.clades(tr2, btr2, rooted=F)


 plot(tr2);title(With sample order randomisation)
 nodelabels(tr2$node.label)

# plot is now ok!! Always 100% for bifurcations and lower for multifurcations!




 --
 Dr. Alastair J. Potts
 Claude Leon Postdoctoral Fellow
 Botany Department
 Nelson Mandela Metropolitan University

 Cell #: 082 491-7275
 Office #: 041 504-2398 (Mon, Wed, Fri)

 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo


___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] R package for detecting positive selection?

[Klaus Schliep]

Dear Anna,

just adding to Emmanuels comment, have a look in the section 4 of the
vignette(phangorn-specials) for a short example.

Cheers,
Klaus


On 6/28/12, Emmanuel Paradis emmanuel.para...@ird.fr wrote:
 Hi Anna,

 See the package phangorn: recent versions implement codon-based models.

 Cheers,

 Emmanuel

 Anna Syme wrote on 28/06/2012 08:42:
 Hi all,
 I was wondering if there are any packages in R that can detect positive
 selection, comparable to programs such as PAML. I have two gene copies
 (each copy in multiple species), and want to see whether either copy has
 been subject to positive selection since divergence.
 Thanks,
 Anna
 *Dr Anna Syme* Molecular Systematist, National Herbarium of Victoria,
 Royal Botanic Gardens Melbourne, Private Bag 2000, Birdwood Avenue,
 South Yarra, 3141, Australia. Phone (03) 9252 2328, Fax (03) 9252 2413.
 anna.s...@rbg.vic.gov.au mailto:anna.s...@rbg.vic.gov.au
 www.rbg.vic.gov.au http://www.rbg.vic.gov.au/

 --
 This email and any attachments may contain information that is personal,
 confidential, legally privileged and/or copyright. No part of it should be
 reproduced, adapted or communicated without the prior written consent of
 the sender and/or copyright owner.

 It is the responsibility of the recipient to check for and remove
 viruses.

 If you have received this email in error, please notify the sender by
 return email, delete it from your system and destroy any copies. You are
 not authorised to use, communicate or rely on the information contained in
 this email.

 Please consider the environment before printing this email.

 This email was Anti Virus checked by RBG Astaro Mail Gateway.



 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

 --
 Emmanuel Paradis
 IRD, Jakarta, Indonesia
 http://ape.mpl.ird.fr/

 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo



-- 
Klaus Schliep
Université Paris 6 (Pierre et Marie Curie)
9, Quai Saint-Bernard, 75005 Paris

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] Cophenetic repeats

[Klaus Schliep]

Hello Ben,

is your tree ultrametric? Do you have a e.g. UPGMA tree?  This would
explain your observation. You can test your tree with
is.ultrametric(trx).

Regards,
Klaus



On 5/27/12, Ben Weinstein bwein...@life.bio.sunysb.edu wrote:
 Hi all,

 I'm trying to decide if this an R error, or an error in
 how I've implemented branch lengths.

 When i create the cophenetic matrix for my tree, and look at the
 relatedness of all tips to a single tip (i.e just looking at one column in
 the matrix). I find that a large portion of them have identical cophenetic
 distances.

  I thought the cophenetic would be the sum of the branch length between
 tips (patristic distance), therefore there should only be pairs (sister
 species have equal terminal branch length)  with *exactly * the same value.
 Comparisons from distantly related taxa should be similar, since most of
 the distance is dictated by the internal branches, but the terminal branch
 should atleast create some difference.

 Sorry that i can't really create a reproducible example for the question.

 Here is some documentation:

 trx-read.nexus(ColombiaPhylogenyUM.tre)
 seeds-sample(trx$tip.label,1)
 sp.weight.alpha-cophenetic(trx)[,seeds]


 table(sp.weight.alpha)sp.weight.alpha
  0 0.00775074  0.1179723  0.2187406  0.2379734 0.23797341
 0.2525792
  1  1  2  4  3  2
 6
 0.36612256 0.36612257 0.52326843 0.59034104 0.59034105 0.59034106
 0.607038
  2  2  2 31 82 29
 3



 82 of the species have the exact same distance to my selected tip .59034105


 Am i using this function correctly?

 I appreciate the help,

 Ben Weinstein

 --
 Ben Weinstein
 Graduate Student
 Ecology and Evolution
 Stony Brook University

 http://life.bio.sunysb.edu/~bweinste/index.html

   [[alternative HTML version deleted]]

 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo


___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] center of tree

[Klaus Schliep]

Hello Mathias,

I'm not aware of any general  COT in R, your chance to make a contribution.
Midpoint rooting is the COT for a special function F [2], it minimizes
the maximum
distance from the root (COT) to the tips.

Cheers,
Klaus



On 5/18/12, Walter, Mathias math...@walter.gd wrote:
 Hi,

 does anyone know a package which computes the center of a tree (COT)?

 The COT is a point on an unrooted phylogeny, where the average
 evolutionary distance from this point to each tip on the phylogeny is
 minimized [1,2].

 Once I thought this COT is identical with the midpoint used to root a
 tree. But I'm unsure.


 [1] http://www.sciencemag.org/content/299/5612/1515.3
 [2] http://jvi.asm.org/content/81/16/8507

 --
 Regards,
 Mathias

 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo



-- 
Klaus Schliep
Université Paris 6 (Pierre et Marie Curie)
9, Quai Saint-Bernard, 75005 Paris

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] treefinder chronogram problems

[Klaus Schliep]

Hello John,

it seems the chronogram contains the heights instead of the
edge.lengths for each edge.
Here is a small function to convert the tree. You should check if it
is the same as in treefinder.

convert.tree - function(tree){
   el = numeric(max(tree$edge))
   el[tree$edge[,2]] = tree$edge.length
   el[min(tree$edge[,1])] = tree$root.edge
   tree$edge.length = el[tree$edge[,1]] - el[tree$edge[,2]]
   tree
}

tree = read.tree(chrono3.tree)
tree2 = convert.tree(tree)
plot(tree2)

Regards,
Klaus


On 4/26/12, John Denton jden...@amnh.org wrote:
 Hi folks,

 I generated a quick chronogram in TreeFinder (attached) for an exploratory
 analysis in R. The tree appeared fine in the Treefinder plot window, but the
 branch length notation does not appear to be compatible with R plotting for
 a chronogram.

 The basal divergence time is 33.9 MYA (which appears in the Newick file),
 but FigTree is telling me the base is something around 22. Has anyone else
 encountered the problem with TreeFinder chronograms? Does anybody have a
 suggestion for how to rescale? I'd prefer to use this tree rather than
 simply make the original unrooted tree ultrametric.

 Thanks!

 ~John


 John S. S. Denton
 Ph.D. Candidate
 Department of Ichthyology and Richard Gilder Graduate School
 American Museum of Natural History
 www.johnssdenton.com

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] parsimony bootstrapping

[Klaus Schliep]

Hello Brian,

it is quite easy:

trees = bootstrap.phyDat(primates, pratchet, trace = 0)
plotBSNew(treeRatchet, trees)

You can use the same parameters as for pratchet and there is some
support for parallel code if you run it on a terminal (on Linux or
Mac).
Unfortunately there is was a bug introduced recently into plotBS, so
here is a working replacement:



plotBSNew - function (tree, BStrees, type = unrooted, bs.col = black,
   bs.adj = NULL, ...)
{

prop.clades - function (phy, ..., part = NULL, rooted = FALSE)
{
  if (is.null(part)) {
  obj - list(...)
  if (length(obj) == 1  class(obj[[1]]) != phylo)
  obj - unlist(obj, recursive = FALSE)
  part - prop.part(obj, check.labels = TRUE)
  }

  bp - prop.part(phy)
  if (!rooted){
  bp - postprocess.prop.part(bp)
  part - postprocess.prop.part(part)   # This line I added!!
  }
  n - numeric(phy$Nnode)
  for (i in seq_along(bp)) {
  for (j in seq_along(part)) {
  if (identical(bp[[i]], part[[j]])) {
  n[i] - attr(part, number)[j]
  done - TRUE
  break
  }
  }
  }
  n
}

   if (type == phylogram | type == cladogram) {
   if (!is.rooted(tree))
   tree2 = midpoint(tree)
   plot(tree2, type = type, ...)
   }
   else plot(tree, type = type, ...)
   x = prop.clades(tree, BStrees)
   x = round((x/length(BStrees)) * 100)
   label = c(rep(, length(tree$tip)), x)
   ind - get(last_plot.phylo, envir = .PlotPhyloEnv)$edge[,
   2]
   if (type == phylogram | type == cladogram) {
   root = getRoot(tree)
   label = c(rep(0, length(tree$tip)), x)
   label[root] = 0
   ind2 = matchEdges(tree2, tree)
   label = label[ind2]
   ind = which(label  0)
   if (is.null(bs.adj))
   bs.adj = c(1, 1)
   nodelabels(text = label[ind], node = ind, frame = none,
   col = bs.col, adj = bs.adj, ...)
   }
   else {
   if (is.null(bs.adj))
   bs.adj = c(0.5, 0.5)
   edgelabels(label[ind], frame = none, col = bs.col,
   adj = bs.adj, ...)
   }
}



Rgeards,
Klaus



On 4/3/12, Brian Bourke brianbou...@hotmail.com wrote:

 Dear All
 I am constructing maximum parsimony trees:

 library(phangorn)
 primates = read.phyDat(primates.dna, format=phylip, type=DNA)
 dm = dist.dna(as.DNAbin(primates))
 treeNJ = NJ(dm)
 treePars = optim.parsimony(treeNJ, primates)
 ##and
 treeRatchet = pratchet(primates, trace = 0)

 , and I would like to know how to do parsimony bootstrap analysis.  Anyone
 know how this can be done?

 Regards,
 Brian.
   
   [[alternative HTML version deleted]]

 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo


___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] Zero bootstrap values

[Klaus Schliep]

Dear Gwennaël,

the bootstrapping function is fine, it is just a problem with the
plotting. Can you try this function, it should work:


plotBSNew - function (tree, BStrees, type = unrooted, bs.col = black,
bs.adj = NULL, ...)
{

prop.clades - function (phy, ..., part = NULL, rooted = FALSE)
{
   if (is.null(part)) {
   obj - list(...)
   if (length(obj) == 1  class(obj[[1]]) != phylo)
   obj - unlist(obj, recursive = FALSE)
   part - prop.part(obj, check.labels = TRUE)
   }

   bp - prop.part(phy)
   if (!rooted){
   bp - postprocess.prop.part(bp)
   part - postprocess.prop.part(part)   # This line I added!!
   }
   n - numeric(phy$Nnode)
   for (i in seq_along(bp)) {
   for (j in seq_along(part)) {
   if (identical(bp[[i]], part[[j]])) {
   n[i] - attr(part, number)[j]
   done - TRUE
   break
   }
   }
   }
   n
}

if (type == phylogram | type == cladogram) {
if (!is.rooted(tree))
tree2 = midpoint(tree)
plot(tree2, type = type, ...)
}
else plot(tree, type = type, ...)
x = prop.clades(tree, BStrees)
x = round((x/length(BStrees)) * 100)
label = c(rep(, length(tree$tip)), x)
ind - get(last_plot.phylo, envir = .PlotPhyloEnv)$edge[,
2]
if (type == phylogram | type == cladogram) {
root = getRoot(tree)
label = c(rep(0, length(tree$tip)), x)
label[root] = 0
ind2 = matchEdges(tree2, tree)
label = label[ind2]
ind = which(label  0)
if (is.null(bs.adj))
bs.adj = c(1, 1)
nodelabels(text = label[ind], node = ind, frame = none,
col = bs.col, adj = bs.adj, ...)
}
else {
if (is.null(bs.adj))
bs.adj = c(0.5, 0.5)
edgelabels(label[ind], frame = none, col = bs.col,
adj = bs.adj, ...)
}
}


plotBSNew(Cassidinaealnx2pmloptim $tree, Cassidinaepmlbs, type=p)


Regards,
Klaus

On 4/3/12, Gwennaël Bataille gwennael.batai...@uclouvain.be wrote:
 Dear all,
 I have got some troubles with the ML bootstrapping {phangorn} : I get
 zero values...
 I already had this problem before, but cannot see what's going wrong here.

 The script and fasta file are in the end of the mail.
 Any idea ?

 Many thanks,


 Gwennaël



 _

 #Import the data  make a first NJ tree
 Cassidinaealnx2- read.dna(file = Cassidinaealnx.fasta, format =
 fasta)
 Cassidinaealnxdist- dist.dna(Cassidinaealnx2, model = K80,
 pairwise.deletion = FALSE)
 Cassidinaetrx- nj(Cassidinaealnxdist)
 plot(Cassidinaetrx)
 Cassidinaertrx- root(Cassidinaetrx, Imatidiumthoracicum)

 #ML
 Cassidinaealnx2pml- pml(Cassidinaertrx, as.phyDat(Cassidinaealnx2))
 Cassidinaealnx2pmloptim- optim.pml(Cassidinaealnx2pml, optNni=TRUE)
 plot(Cassidinaealnx2pmloptim)
 Cassidinaepmlbs- bootstrap.pml(Cassidinaealnx2pmloptim, optNni=TRUE)
 plotBS(Cassidinaealnx2pmloptim $tree, Cassidinaepmlbs, type=p)




 
 Cassidinaealnx.fasta





Agroiconotapropinqua
 tganaaaccgaa-ggatcg-aataaattcattcgcg-tttcgattgtcagcgttgagcgg-ttgtgtgactgacggagtacgcttcg-tccgcaccttcagttcatcga-act-ggcttcgaacgcgtgcac-tagtaggacgtcgcgatccg-ttgggcgtcggtctaaggcttgcgg-tggagctcgcgtg--gacgtttc--ggcgtt--ccgcacgaacccgcaatgtccc-gatcgactcgctcgacggt-ata-aagatggcgcccgctactcagttagcgttcgacctacggca-agc-acgtccggtgttt--gatggcn-atcagacctggtcggatcctgtcctc-ggacgactgttgg--ctcgtantg-ttctcgaacagacctcg-ttkaaacgcc-catctgcgacgcgctt-tgggtactttca-ggacccgtcttgaaacacgg-
Microctenochiracumulata
 ---ctgagaamccgaa-agatcg-aataaattcattcgcg-tttcgattgtcagcgttgagcgg-ttgtgtgactgacggagtacgcttcg-tccgcaccttcagttcatcga-act-ggcttcgaacgcgtgcac-tagtaggacgtcgcgatccg-ttgggcgtcggtctaaggcttgcgg-tggagctcgcgtg--gacgtttc--ggcgt---ccgcacgaacccgcaatgtccc-gatcgacacgctcgacggt-aca-aagatggcgcccgctattcagttagcgttcgacctatggca-agc-acgtccggtgttt--gatggcg-atcagacctggtcggatcccgtcctc-ggacgactgttgg--ctcgtagtg-ttctcgaacagacctcg-ttgaaacgcc-gatctgcgacgctatagctt-tgggtactttca-ggacccgtcttgaaacacg--
Charidotellasexpunctata
 --a-tcg-aataaattcatkggcg-tttcgattgtcagcgttggacgg-ttgtgtgactgacggagtacgcttcg-tccgcaccttcagttcatcga-act-ggcttcgaacgcgtgcac-tagtaggacgtcgcgatccg-ttgggcgtcgatctaaggcttgcgg-tggagctcgcgtg--gacgtttc--ggcgtt--ccacacgaacctgcaatgtccc-gatcgactcgctcgacggt-ata-aagatggcgcccgctattcagttagcgttcggcctacggca-agc-acgtgcggtattt--gatggcg-atcagacctggtcggatcccgtcgcc-ggacgactgttgg--ctcgtagtg-ttctcgaacagacctcg-ttgaaacgcc-gatctgcgacgctntagctt-tgggtactttca-ggacccgtcttgaaacacggcagg-
Charidotellasinuata
 

Re: [R-sig-phylo] Boostrap and likelihood values on branches

[Klaus Schliep]

Dear Gwennaël,

first it seems there was a bug in the function prop.clades from ape
introduced recently, results in many zeros of bootstrap values.
This function needs to be replaced in the ape package:

prop.clades - function (phy, ..., part = NULL, rooted = FALSE)
{
if (is.null(part)) {
obj - list(...)
if (length(obj) == 1  class(obj[[1]]) != phylo)
obj - unlist(obj, recursive = FALSE)
part - prop.part(obj, check.labels = TRUE)
}

bp - prop.part(phy)
if (!rooted){
bp - postprocess.prop.part(bp)
part - postprocess.prop.part(part)   # This line I added!!
}
n - numeric(phy$Nnode)
for (i in seq_along(bp)) {
for (j in seq_along(part)) {
if (identical(bp[[i]], part[[j]])) {
n[i] - attr(part, number)[j]
done - TRUE
break
}
}
}
n
}



On 3/12/12, Gwennaël Bataille gwennael.batai...@uclouvain.be wrote:
 Dear all,
 I used a very raw 18S phylogeny of Craniates to introduce phylogeny to
 students, and I would like to estimate bootstrap or likelihood values
 for branches. Sorry for this long message, but I hope you can answer at
 least some of my questions.

 1) For bootstrapping, I obtain zero values for some of them and wonder
 : i) if it is possible (and how to interpret it) ; ii) what to change to
 correct it ?
 Here is the script (the alnx.fasta file can be found at the end of
 this message) :

 library(ape)
 alnx2 - read.dna(file = alnx.fasta, format = fasta)
 alnxdist - dist.dna(alnx2, model = raw, pairwise.deletion = FALSE)
 trx - nj(alnxdist)
 rtrx - root(trx, Branchiostoma)
 plot(rtrx)

 bp - boot.phylo(rtrx, alnx2, function(x) nj(dist.dna(x, model = raw,
 pairwise.deletion = FALSE))) # I also obtain zero values if
 pairwise.deletion = TRUE
 nodelabels(bp)


 2) Secondly, I also would like to obtain likelihood values for branches
 (I think it is the second common kind of values found on branches, with
 bootstrap ones), but could not find a function to do so. Would you know
 one ?
 If I try maximum likelihood estimation, are those likelihood values
 calculated for the branches ? I see them for trees and sites, but not
 for branches…

The branch length in a ML tree are (proportional to) time (molecular
clock) or the expected number of substitutions per site.


 library(phangorn)
 alnx2pml - pml(rtrx, as.phyDat(alnx2))
 alnx2pmloptim - optim.pml(alnx2pml) # This is just a test ; I don't
 really need a model
 plot(alnx2pmloptim)
 #The following commands do not seem to be a great idea (for boostrap) ;
 my computer doesn't like them…
 # bp2 - boot.phylo(rtrx, alnx2, function(x) optim.pml( pml(
 nj(dist.dna(x)) , as.phyDat(x)) ) )
 # nodelabels(bp2)


There is function for bootstrapping in phangorn:

library(phangorn)
alnx2pml - pml(rtrx, as.phyDat(alnx2))
alnx2pmloptim - optim.pml(alnx2pml, TRUE)
# some tree rearrangements, otherwise you may get some zeros
alnx2bs - bootstrap.pml(alnx2pmloptim, optNni=TRUE)
plotBS(alnx2pmloptim$tree, alnx2bs)



 3) Then, where to find likelihood values after a mrbayes run ?

 library(phyloch)
 mrbayes(alnx2, path = getwd(),ngen = 100, file = test.con, run =
 FALSE) # I obtain the file I can use then in mrbayes


 4) Eventually, for parsimony, what do the parsimony score represent ? Is
 it a way to weight it ? I suppose the answer is no, and parsimony scores
 are only there to compare them, not to deduce what is a good or bad
 score…

The parsimony score is the minimal number of substitutions needed to
account for the data on a phylogeny.
If you weigh a tree with the parsimony score the edge length represent
the minimal number of substitutions needed.


 alnx2pars - parsimony(rtrx, as.phyDat(alnx2))
 alnx2parsoptim - optim.parsimony(rtrx, as.phyDat(alnx2))
 plot(alnx2parsoptim)

pratchet is likely to find better tree (not for this small tree here)

alnx2pars - parsimony(rtrx, as.phyDat(alnx2))
alnx2parsoptim - pratchet(as.phyDat(alnx2), start=rtrx)
alnx2parsoptim - acctran(alnx2parsoptim, as.phyDat(alnx2))
  # assign edge length via ACCTRAN criterion,
  # the assignment for the edge length is in general not unique!!!

tree - midpoint(alnx2parsoptim) # midpoint rooting for a nicer look
plot(tree)
edgelabels(round(tree$edge.length, 3), adj = c(.5, 1), frame=none)
add.scale.bar()






 Thanks a lot in advance for your answer.
 All the best,


 Gwennaël Bataille

 ___


  Salmo
 

Re: [R-sig-phylo] Estimating branch lengths on a fixed topology from molecular data

[Klaus Schliep]

Hello Alejandro,

you can use the phangorn package to do this in R.
The code could look like this:

library(phangorn)
tree - read.tree(treefile) # load the tree / topology
dat - read.phyDat(datafile) # load the data

#Parsimony, apply length using the ACCTRAN criterion
treeMP - acctran(tree, data)
plot(treeMP)

# ML
fit - pml(tree, data) # set up a ML object
fit - optim.pml(fit) # optimise edge lengths
treeML - fit$tree
plot(treeML)

Regards,
Klaus


On 2/3/12, Alejandro Gonzalez V alejandro.gonza...@ebd.csic.es wrote:
 Hello,

 I was wondering if it is possible to estimate branch lengths on a fixed
 topology from a matrix of nuclear or mitochondrial DNA data. I used to do
 this in PAUP and would very much like to have an alternative in R. In PAUP
 one loaded a matrix of molecular data and a topology (without branch
 lengths) with identical species coverage, then branch lengths were estimated
 using maximum likelihood (could also use parsimony) and given settings for
 the substitution model.
 Any tips would be appreciated!

 Best wishes

 Alejandro
 __

 Alejandro Gonzalez Voyer
 Post-doc

 Estación Biológica de Doñana (CSIC)
 Avenida Américo Vespucio s/n
 41092 Sevilla
 Spain

 Tel: +34- 954 466700, ext 1749

 E-mail: alejandro.gonza...@ebd.csic.es

 Web-site (Under Construction):
 Group page: http://consevol.org/index.html
 Personal web-page: http://consevol.org/members/alejandro.html















   [[alternative HTML version deleted]]



___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] phyDat format problem

[Klaus Schliep]

Hi Jenny,

thanks Liam. he's right. You have to change the data to a matrix.

In a data.frame all elements of a columns must contain the same data
type. When using data.frames I stored characters in the columns. This
would make it possible to allow to have binary characters in one
column and nucleotides or amino acid in others.

However alignments often come in form that in each that the taxa are
stored in rows and characters in columns and phyDat assumes that a
matrix is stored in this form.

The transpose function t(p) has the side effect that it changes a
data.frame into a matrix, that's why your try didn't work out.
 class(p)
[1] data.frame
 class(ptrans)
[1] matrix

Regards,
Klaus



On 1/9/12, Liam J. Revell liamjrev...@gmail.com wrote:
 Hi Jen.

 Try changing to a matrix first (from a data frame). [I.e.,
 p-as.matrix(p); p2-phyDat(p,type=USER,levels=c(0,1)) ]

 I'm not sure why this works, but it seems to. - Liam

 More details below.

 When I do:
 p-matrix(c(1,0,0,0,1,0,1,0,1,0,
 1,1,0,1,1,0,0,1,0,0,
 1,0,0,1,0,1,0,0,1,0,
 0,1,0,1,1,1,0,0,0,0,
 0,0,0,1,1,1,1,0,1,0),5,10,byrow=T)
 dimnames(p)-list(c(A,B,C,D,E),paste(X,1:10,sep=))
 p2-phyDat(p,type=USER,levels=c(0,1))

 I get:
   p2
 5 sequences with 10 character and 9 different site patterns.
 The states are 0 1

 However:
 p-as.data.frame(p)
 p2-phyDat(p,type=USER,levels=c(0,1))

 Gives me:
   p2
 10 sequences with 5 character and 5 different site patterns.
 The states are 0 1

 --
 Liam J. Revell
 University of Massachusetts Boston
 web: http://faculty.umb.edu/liam.revell/
 email: liam.rev...@umb.edu
 blog: http://phytools.blogspot.com


 On 1/9/2012 9:14 AM, J Greenwood wrote:
 Hi all,

 I am having a problem with getting parsimony scores in Phangorn and have
 found that the program is not reading in my datafile the way I expect it.
 I have 10 characters and 5 taxa, but phangorn seems to read this the
 opposite way round.

 Can anybody help? Details follow,

 thanks,

 Jen


 I made a test dataset to try and work out the problem which is the
 following:

 p-read.table(file=tab_phangorn.txt, header=T)
 p
 X1 X2 X3 X4 X5 X6 X7 X8 X9 X10
 A 1 0 0 0 1 0 1 0 1 0
 B 1 1 0 1 1 0 0 1 0 0
 C 1 0 0 1 0 1 0 0 1 0
 D 0 1 0 1 1 1 0 0 0 0
 E 0 0 0 1 1 1 1 0 1 0

 Where A-E are taxa and X1-X10 are characters.
 I then tried converting it to phyDat format which gives the following:

 p2-phyDat(p, type=USER, levels=c(0,1))
 p2

 10 sequences with 5 character and 5 different site patterns.
 The states are 0 1

 str(p2)
 List of 10
 $ X1 : int [1:5] 2 2 2 1 1
 $ X2 : int [1:5] 1 2 1 2 1
 $ X3 : int [1:5] 1 1 1 1 1
 $ X4 : int [1:5] 1 2 2 2 2
 $ X5 : int [1:5] 2 2 1 2 2
 $ X6 : int [1:5] 1 1 2 2 2
 $ X7 : int [1:5] 2 1 1 1 2
 $ X8 : int [1:5] 1 2 1 1 1
 $ X9 : int [1:5] 2 1 2 1 2
 $ X10: int [1:5] 1 1 1 1 1
 - attr(*, class)= chr phyDat
 - attr(*, weight)= int [1:5] 1 1 1 1 1
 - attr(*, nr)= int 5
 - attr(*, nc)= int 2
 - attr(*, index)= int [1:5] 1 2 3 4 5
 - attr(*, levels)= num [1:2] 0 1
 - attr(*, allLevels)= chr [1:3] 0 1 ?
 - attr(*, type)= chr USER
 - attr(*, contrast)= num [1:3, 1:2] 1 0 1 0 1 1

 It has read the data in the opposite way round to the way i was
 expecting, which was odd, but I transposed the matrix anyway, but got
 the same output:

 ptrans-t(p)
 ptrans
 A B C D E
 X1 1 1 1 0 0
 X2 0 1 0 1 0
 X3 0 0 0 0 0
 X4 0 1 1 1 1
 X5 1 1 0 1 1
 X6 0 0 1 1 1
 X7 1 0 0 0 1
 X8 0 1 0 0 0
 X9 1 0 1 0 1
 X10 0 0 0 0 0

 ptrans2
 10 sequences with 5 character and 5 different site patterns.
 The states are 0 1
 str(ptrans2)
 List of 10
 $ X1 : int [1:5] 2 2 2 1 1
 $ X2 : int [1:5] 1 2 1 2 1
 $ X3 : int [1:5] 1 1 1 1 1
 $ X4 : int [1:5] 1 2 2 2 2
 $ X5 : int [1:5] 2 2 1 2 2
 $ X6 : int [1:5] 1 1 2 2 2
 $ X7 : int [1:5] 2 1 1 1 2
 $ X8 : int [1:5] 1 2 1 1 1
 $ X9 : int [1:5] 2 1 2 1 2
 $ X10: int [1:5] 1 1 1 1 1
 - attr(*, class)= chr phyDat
 - attr(*, weight)= int [1:5] 1 1 1 1 1
 - attr(*, nr)= int 5
 - attr(*, nc)= int 2
 - attr(*, index)= int [1:5] 1 2 3 4 5
 - attr(*, levels)= num [1:2] 0 1
 - attr(*, allLevels)= chr [1:3] 0 1 ?
 - attr(*, type)= chr USER
 - attr(*, contrast)= num [1:3, 1:2] 1 0 1 0 1 1

 Can anyone explain to me what I am doing wrong? How can I get the phyDat
 format to hold my information correctly?

 --
 J Greenwood
 jenny.greenw...@bristol.ac.uk

 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo


 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo



-- 
Klaus Schliep
Université Paris 6 (Pierre et Marie Curie)
9, Quai Saint-Bernard, 75005 Paris

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] Error message using clustal{ape} under Windows

[Klaus Schliep]

Dear Gwennaël,

it seems it is a problem with the path to clustal, as the path
contains spaces. Try
clustal(x, exec = shortPathName(C:/Program Files
(x86)/ClustalW2/clustalw2.exe))
it works for me.

Regards,
Klaus



On 1/3/12, Gwennaël Bataille gwennael.batai...@uclouvain.be wrote:
 Dear all,
 I really need some help for using the clustal function of the ape
 package under Windows... I already discussed this with Emmanuel Paradis
 (in charge of this function), but the data seem to work correctly under
 his Linux system, and he was unfortunately not able to fix this problem.

 Here is my basic script :
 library(ape)
 x - read.dna(C:\\Users\\ gwbataille\\Documents\\ ...
 \\Test_Plusieurs_Sequences_EmAlbicoxaget_New_Method.fasta, format =
 fasta)
 clustal(x)

 And the error message is :
 Error in file(file, r) : cannot open the connection
 In addition: Warning message:
 In file(file, r) :
 cannot open file
 'C:\Users\GWBATA~1\AppData\Local\Temp\RtmpE3v3p7/input_clustal.aln':
 No such file or directory

 I also tried :
 clustal(x, exec = C:/Program Files (x86)/ClustalW2/clustalw2.exe)

 and :
 fix(clustal)
 ## change outf - paste(d, input_clustal.aln, sep = /) - sep = \\

 But always get the same error message…


 Thank you very much in advance for your help.

 Regards,


 Gwennaël Bataille

 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo


___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] basal dichotomy in phylo object

[Klaus Schliep]

Dear Ondrej,

there is also a function midpoint in phangorn, which roots the trees
in the middle of the longest path. Maybe you are lucky and your edge
(500, 174) is one in the middle.

Regards,
Klaus

On 10/25/11, Liam J. Revell liam.rev...@umb.edu wrote:
 In my phytools package
 (http://cran.r-project.org/web/packages/phytools/index.html) there is a
 function reroot() that is really just a simple wrapper around root()
 that allows you to root along internal or terminal edges (rather than
 only at nodes).

 In your example, assuming that the node 500 is tipward of 174 in the
 tree, to re-root the tree halfway along the edge leading to node number
 500 (in tree$edge), you would just do:

 rooted.tree-reroot(tree,node.number=500,position=tree$edge.length[tree$edge[,2]==500]/2)

 I think.

 All the best, Liam

 --
 Liam J. Revell
 University of Massachusetts Boston
 web: http://faculty.umb.edu/liam.revell/
 email: liam.rev...@umb.edu
 blog: http://phytools.blogspot.com

 On 10/25/2011 10:44 AM, Ondřej Mikula wrote:
 Hello,
 I have an unrooted tree, but I know that basal dichotomy is between nodes
 labeled 174 and 500. I am wandering how to indicate it in the object
 of
 class 'phylo'. I tried 'root' function of 'ape' but with no success.
 I will be grateful for any advice. Best wishes
 Ondrej Mikula


 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo



-- 
Klaus Schliep
Université Paris 6 (Pierre et Marie Curie)
9, Quai Saint-Bernard, 75005 Paris

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] Negative Eigenvalues from Phylo Distance Matrix

[Klaus Schliep]

Hi Jonathan,

you probably mixed up positive matrices and positive definite matrices.
To have only positive eigen values a matrix must be positive definite.
Covariance matrices are for example are always positive semi-definite,
even though they can contain negative entries.
D is symmetric therefore all eigen values are reel and positive
therefore the largest eigen value is positive (Perron–Frobenius
theorem).

Regards
Klaus

On 9/27/11, Jonathan Mitchell jonsmit...@gmail.com wrote:
 Hello!

 I tried on some simulated trees last night and they weren't working, but
 after your message I started a fresh R session and now solve() works fine
 for both my ultrametric tree, and my simulated trees. No idea what happened
 there, but thanks!

 Unfortunately, I still am getting almost all negative eigenvalues for the
 decomposition of the distance matrix, which doesn't make sense to me. I was
 under the impression that the distance matrix from an ultrametric tree
 should be Euclidean. Further, I can confirm that the distance matrix is
 supposed to be euclidean by using the is.euclid() function from ade4
 [is.euclid(as.dist(D)) returns TRUE], which means I have no idea why
 eigen(D) is giving me negative eigenvalues...

 Any thoughts?

 --Jon

 On Tue, Sep 27, 2011 at 12:23 AM, Emmanuel Paradis
 emmanuel.para...@ird.frwrote:

 Hi John,

 I tried with a couple of simulated trees (ultrametric or not) and got no
 error. Have you tried using the 'tol' argument of solve()?

 Emmanuel
 -Original Message-
 From: Jonathan Mitchell mitchel...@uchicago.edu
 Sender: r-sig-phylo-boun...@r-project.org
 Date: Mon, 26 Sep 2011 17:42:49
 To: r-sig-phylo@r-project.org
 Reply-To: jonsmit...@gmail.com
 Subject: [R-sig-phylo] Negative Eigenvalues from Phylo Distance Matrix

 Hello all,

 I've run into a problem trying to find the eigen vectors for a
 phylogenetic
 distance matrix. Namely, the transposed vector matrix crossed with the
 matrix itself is singular. Also, the eigenvalues for the distance matrix
 proper are almost all negative, which I don't understand whatsoever.

 #Given a tree with branch lengths (in my case, it's ultrametric)
 D - cophenetic(tree)
 eD - eigen(D, symmetric=TRUE)$vectors  #note that all elements except the
 first in eD$values are negative for me
 inv - solve(t(eD)%*%eD)  #this is what I really want, the inverse of the
 transposed matrix of vectors multiplied by itself

 Error in solve.default(t(eD) %*% eD) :
  system is computationally singular: reciprocal condition number =
 2.46807e-29

 Any help would be greatly appreciated!

 -- Jon

 _

 Jonathan S. Mitchell
 http://home.uchicago.edu/~mitchelljs/

 PhD Student
 Committee on Evolutionary Biology
 The University of Chicago
 1025 57th Str, Culver Hall 402
 Chicago, IL 60637

 Geology Department
 The Field Museum of Natural History
 1400 S. Lake Shore Dr.
 Chicago, IL 60605

 [[alternative HTML version deleted]]

 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo


   [[alternative HTML version deleted]]

 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo



-- 
Klaus Schliep
Université Paris 6 (Pierre et Marie Curie)
9, Quai Saint-Bernard, 75005 Paris

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] phylogenetic tree diameter

[Klaus Schliep]

Hi John,

On 7/28/11, John Denton jden...@amnh.org wrote:
 Hi all,

 I'm using the ape package, and was wondering if there was a method in the
 package or code for calculating the diameter

 (http://mathworld.wolfram.com/GraphDiameter.html)

 of phylogenetic tree objects. I had been looking in to exporting trees
 from ape into a graph theory package (igraph), but it looks like it would
 be more difficult to write a file conversion script than to try something
 in ape itself.

Conversion is not hard, look in the source code of
phangorn:::plot.networx  for an example.


 Also, I have tried installing the phybase package from CRAN but keep
 getting errors. Does anyone have or know of a good script for calculating
 Robinson-Foulds differences for phy objects?

dist.topo in ape and RF.dist and treedist in phangorn compute the
Robinson-Foulds or symmetric distance (called Penny-Hendy in
dist.topo).

Cheers,
Klaus

 Thanks!

 ~J



 --
 ORGANIC LIFE beneath the shoreless waves
 Was born and nurs'd in ocean's pearly caves;
 First forms minute, unseen by spheric glass,
 Move on the mud, or pierce the watery mass;
 These, as successive generations bloom,
 New powers acquire and larger limbs assume;
 Whence countless groups of vegetation spring,
 And breathing realms of fin and feet and wing.

 ~Erasmus Darwin, The Temple of Nature Canto I.V

 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo



-- 
Klaus Schliep
Université Paris 6 (Pierre et Marie Curie)
9, Quai Saint-Bernard, 75005 Paris

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] phylo to dendrogram

[Klaus Schliep]

Hello Jarrod,

try:
h.d-as.hclust(reorder(multi2di(host.tree)))
as.dendrogram(h.d)
Maybe it works.

Regards,
Klaus


On 7/25/11, Jarrod Hadfield j.hadfi...@ed.ac.uk wrote:
 Dear list,

 I'm trying to convert a phylo object into a dendrogram object with
 little success. host.tree is a rooted ultrametric tree with polytomies
 stored as phylo object. The polytomies can be resolved to produce a
 binary tree and this can be coerced into a hclust object:

 h.d-as.hclust(multi2di(host.tree))

 However when converted to a dendrogram using as.dendrogram(h.d) the
 dendrogram appears broken. str(as.dendrogram(h.d)) returns:

 --[dendrogram w/ 2 branches and  members at h = 96.1]
|--[dendrogram w/ 2 branches and  members at h = 82.8]
|  |--[dendrogram w/ 2 branches and  members at h = 70.3]
|  |  |--[dendrogram w/ 2 branches and  members at h = 39]
|  |  |  |--[dendrogram w/ 2 branches and  members at h = 24.5]
 Error in class(val) - cl : attempt to set an attribute on NULL

 I have given the zero branch lengths from the polytomies non-zero
 branch lengths, but this does not seem to work either.

 Any help would be gratefully received.

 Jarrod Hadfield


 --
 The University of Edinburgh is a charitable body, registered in
 Scotland, with registration number SC005336.

 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo



-- 
Klaus Schliep
Université Paris 6 (Pierre et Marie Curie)
9, Quai Saint-Bernard, 75005 Paris

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] outputting values during taxon pruning

[Klaus Schliep]

Hi John

I added a few comments.

Regards,
Klaus

On 7/25/11, John Denton jden...@amnh.org wrote:
 Hi, all.

 I'm working on a script that iteratively removes taxa from an initial
 total tree and re-searches for optimal (pruned) topologies after
 outputting a pruned alignment.

 I'm trying to get a few things to happen:

 First, the script iteratively checks to make sure that the taxon names
 match in the alignments and trees, and terminates if they do not. I can
 get it to print that the matches are okay as it checks each taxon:

 for (i in totaltr$tip.label) {
   if (all(totaltr$tip.label %in% rownames(alig))  all(rownames(alig) 
 %in%
 totaltr$tip.label)) cat(OK: Taxon labels match in tree and alignment)
   else stop(Taxon labels do not match. Script terminates.)
   }

 But I want each line of OK: Taxon labels match and Taxon labels do
 not match to print the name of the taxon being checked, or the taxon for
 which a match does not occur, something like

 Taxon labels match in tree and alignment FOR TAXON i. and

 Taxon labels do not match FOR TAXON i. Script terminates.

 Suggestions?

use paste instead of cat, you have done it below.



 Second, I want the script to store a list of the tree lengths it
 calculates as it goes. I can calculate tree lengths via
 sum(tree.name$edge.length), but I am having trouble, it seems,
 initializing and storing the list. The code I have is

 ts.values - list()
 length(ts.values) - length(row.names(alig))

that doesn't work, use
ts.values - vector(list, length(row.names(alig)))
but a numeric/double vector would probably easier to handle




 for (i in totaltr$tip.label) {
   if (all(totaltr$tip.label %in% rownames(alig))  all(rownames(alig) 
 %in%
 totaltr$tip.label)) cat(OK: Taxon labels match in tree and alignment)
   else stop(Taxon labels do not match. Script terminates.)
   }

 Sys.sleep(3)

k=1
 for (i in totaltr$tip.label) {
# i is now a tip.label
   write.tree(drop.tip(totaltr, i), file = paste(minus_, i, .tre, sep =
 ) )
   
# here you need a number or you have to set names before
   ts.values[[k]] -  sum(drop.tip(totaltr, i)$edge.length)
k=k+1
 }

You can use ts.values[[i]] if you set names of beforehand
attr(ts.values, names) = row.names(alig)



 So ts.values is the list I am trying to store the tree lengths in. After
 the script finishes, though, I try to call ts.values via print(ts.values),
 but it says the object does not exist.

 I am tinkering with these two problems using a modified version of the
 phymltest routine in the ape package, so the tree/alignment filenames
 written are for phyml output.

 Suggestions?

Have you used phangorn? It also contains a function modelTest, so you
may give it a try.


 Thanks a lot.

 ~John




 --
 ORGANIC LIFE beneath the shoreless waves
 Was born and nurs'd in ocean's pearly caves;
 First forms minute, unseen by spheric glass,
 Move on the mud, or pierce the watery mass;
 These, as successive generations bloom,
 New powers acquire and larger limbs assume;
 Whence countless groups of vegetation spring,
 And breathing realms of fin and feet and wing.

 ~Erasmus Darwin, The Temple of Nature Canto I.V

 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo



-- 
Klaus Schliep
Université Paris 6 (Pierre et Marie Curie)
9, Quai Saint-Bernard, 75005 Paris

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] Dating trees using date.phylo - ape problem?

[Klaus Schliep]

Dear Roland,

it would be good if you add the datasets, that one can reproduce your
results (archotreeresolved, ages). I would guess from the error
message that your ages list may has the wrong format.
Also traceback() sometimes tells you where the error happens exactly.

Kind regards
Klaus

On 7/21/11, Roland Sookias r.sook...@gmail.com wrote:
 Hi all

 I'm trying to assign lengths to branches in a tree using the Ruta et al.
 2008 method via Graeme Lloyd's date.phylo() function (see
 http://www.graemetlloyd.com/methdpf.html). I keep coming up against an error
 for some reason. Graeme himself does exactly the same thing seemingly on his
 Mac (I'm using Windows 7, R 2.13.1, Ape 2.7-2 - Graeme has Ape 2.7-1) and it
 works. I've tried various permutations of the tree, with and without a
 rooted symbol [r], and there are no special characters etc., the names in
 tree and ages list are the same. I still keep getting:

  ttree-date.phylo(archotreeresolved, ages, rlen=1, method=equal)
 Error in ages[tree$tip.label[find.descendants(nodes[i], tree)], 1] :
   incorrect number of dimensions

 Any help, or advice what to try/look into much appreciated. Graeme thinks
 it's something with ape not date.phylo...

 Very best

 Roland

   [[alternative HTML version deleted]]

 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo



-- 
Klaus Schliep
Université Paris 6 (Pierre et Marie Curie)
9, Quai Saint-Bernard, 75005 Paris

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] Dating trees using date.phylo - ape problem?

[Klaus Schliep]

Hi Roland,

as I suspected, your ages list was not in the right format (has
nothing to do with ape).
Try this:

tmp = read.csv(agescut.csv, header=FALSE)
ages = matrix(tmp[, 2], ncol=1)
rownames(ages) = tmp[,1]
tree = read.tree(archotreeresolved6.tre)
ttree-date.phylo(tree, ages, rlen=1, method=equal)


Regards,
Klaus



On 7/21/11, Roland Sookias r.sook...@gmail.com wrote:
 Hi

 Thanks guys.

 Attached are the two files (will these work via the list?)

 When I type trackback() all I get is 1: date.phylo(archotreeresolved, ages,
 rlen = 1, method = equal).

 Roland



 On Thu, Jul 21, 2011 at 1:24 PM, Klaus Schliep
 klaus.schl...@gmail.comwrote:

 Dear Roland,

 it would be good if you add the datasets, that one can reproduce your
 results (archotreeresolved, ages). I would guess from the error
 message that your ages list may has the wrong format.
 Also traceback() sometimes tells you where the error happens exactly.

 Kind regards
 Klaus

 On 7/21/11, Roland Sookias r.sook...@gmail.com wrote:
  Hi all
 
  I'm trying to assign lengths to branches in a tree using the Ruta et al.
  2008 method via Graeme Lloyd's date.phylo() function (see
  http://www.graemetlloyd.com/methdpf.html). I keep coming up against an
 error
  for some reason. Graeme himself does exactly the same thing seemingly on
 his
  Mac (I'm using Windows 7, R 2.13.1, Ape 2.7-2 - Graeme has Ape 2.7-1)
  and
 it
  works. I've tried various permutations of the tree, with and without a
  rooted symbol [r], and there are no special characters etc., the names
 in
  tree and ages list are the same. I still keep getting:
 
   ttree-date.phylo(archotreeresolved, ages, rlen=1, method=equal)
  Error in ages[tree$tip.label[find.descendants(nodes[i], tree)], 1] :
incorrect number of dimensions
 
  Any help, or advice what to try/look into much appreciated. Graeme
  thinks
  it's something with ape not date.phylo...
 
  Very best
 
  Roland
 
[[alternative HTML version deleted]]
 
  ___
  R-sig-phylo mailing list
  R-sig-phylo@r-project.org
  https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
 


 --
 Klaus Schliep
 Université Paris 6 (Pierre et Marie Curie)
 9, Quai Saint-Bernard, 75005 Paris




-- 
Klaus Schliep
Université Paris 6 (Pierre et Marie Curie)
9, Quai Saint-Bernard, 75005 Paris

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] Dating trees using date.phylo - ape problem?

[Klaus Schliep]

Hi Roland!
ages should be a matrix, with row names and the ages in the first column.

On 7/21/11, Roland Sookias r.sook...@gmail.com wrote:
 Great! It worked. What format should the file be in then? I thought I'd
 followed the format specified.. Thanks very much indeed :)

 On Thu, Jul 21, 2011 at 4:10 PM, Klaus Schliep
 klaus.schl...@gmail.comwrote:

 Hi Roland,

 as I suspected, your ages list was not in the right format (has
 nothing to do with ape).
 Try this:

 tmp = read.csv(agescut.csv, header=FALSE)
 ages = matrix(tmp[, 2], ncol=1)
 rownames(ages) = tmp[,1]
 tree = read.tree(archotreeresolved6.tre)
 ttree-date.phylo(tree, ages, rlen=1, method=equal)


 Regards,
 Klaus



 On 7/21/11, Roland Sookias r.sook...@gmail.com wrote:
  Hi
 
  Thanks guys.
 
  Attached are the two files (will these work via the list?)
 
  When I type trackback() all I get is 1: date.phylo(archotreeresolved,
 ages,
  rlen = 1, method = equal).
 
  Roland
 
 
 
  On Thu, Jul 21, 2011 at 1:24 PM, Klaus Schliep
  klaus.schl...@gmail.comwrote:
 
  Dear Roland,
 
  it would be good if you add the datasets, that one can reproduce your
  results (archotreeresolved, ages). I would guess from the error
  message that your ages list may has the wrong format.
  Also traceback() sometimes tells you where the error happens exactly.
 
  Kind regards
  Klaus
 
  On 7/21/11, Roland Sookias r.sook...@gmail.com wrote:
   Hi all
  
   I'm trying to assign lengths to branches in a tree using the Ruta et
 al.
   2008 method via Graeme Lloyd's date.phylo() function (see
   http://www.graemetlloyd.com/methdpf.html). I keep coming up against
 an
  error
   for some reason. Graeme himself does exactly the same thing seemingly
 on
  his
   Mac (I'm using Windows 7, R 2.13.1, Ape 2.7-2 - Graeme has Ape 2.7-1)
   and
  it
   works. I've tried various permutations of the tree, with and without
   a
   rooted symbol [r], and there are no special characters etc., the
 names
  in
   tree and ages list are the same. I still keep getting:
  
ttree-date.phylo(archotreeresolved, ages, rlen=1, method=equal)
   Error in ages[tree$tip.label[find.descendants(nodes[i], tree)], 1] :
 incorrect number of dimensions
  
   Any help, or advice what to try/look into much appreciated. Graeme
   thinks
   it's something with ape not date.phylo...
  
   Very best
  
   Roland
  
 [[alternative HTML version deleted]]
  
   ___
   R-sig-phylo mailing list
   R-sig-phylo@r-project.org
   https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
  
 
 
  --
  Klaus Schliep
  Université Paris 6 (Pierre et Marie Curie)
  9, Quai Saint-Bernard, 75005 Paris
 
 


 --
 Klaus Schliep
 Université Paris 6 (Pierre et Marie Curie)
 9, Quai Saint-Bernard, 75005 Paris




-- 
Klaus Schliep
Université Paris 6 (Pierre et Marie Curie)
9, Quai Saint-Bernard, 75005 Paris

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] Bootstrap values and NJ when there is no genetic distance between samples

[Klaus Schliep]

Hi Alastair,

it is not that surprising. NJ normally does not produce poytomies,
just edge weights of length 0. How these are broken may depends from
the input order (from labels in the distance matrix like in this
implementation) or could be broken randomly.  I added some code below
to highlight it.
However you can use the function di2multi to produce real
multifurcating trees.

Cheers,
Klaus

On 5/5/11, Alastair Potts pott...@gmail.com wrote:
 Good day all,
 I noticed something that I would consider an anomaly when analysing one
 of my trees with NJ.

 A 'polytomy' of samples contained many bootstrap values of 100 between
 samples. I was looking at the total change in bootstrap values for all
 nodes when I picked this up (as the signal in the data started to
 disappear during subsetting, the number of nodes with bootstrap values
 of 100 suddenly rocketed up).

 Here's my best shot at reproducing my results in a standard way.

 library(ape)
 a - as.DNAbin(matrix('a',10,10)) # DNA data with no variation
 b - dist.dna(a,mode='N')
 c - nj(b)
 plot(c) #No surprises - a single polytomy

That is wrong the tree has all internal edges, just the edge length are 0.
If you count the tree is fully resolved
 c$edge.length
 [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
c$edge.length[] = 1 # I don't add edges, just change the edge length
plot(c)

 boot.phylo(tr, a, function(xx)
 nj(dist.dna(xx[sample(1:nrow(xx)),],model=N)),B=100)
 # Not what I expected - I would have thought all BS values would be 0

 Am I missing something about how NJ is working in ape. If I perform a
 similar analysis in PAUP none of the bipartions are greater than 12%,
 which is what I would have expected, considering that there is no
 information in the dataset.

 Any help would be greatly appreciated.
 Cheers,
 Alastair




 --
 -
 Alastair Potts
 PhD candidate
 Botany Department
 University of Cape Town
 alastair.po...@uct.ac.za or pott...@gmail.com
 University Private Bag, Rondebosch 7700, South Africa
 or
 PO Box 115, Loxton 6985, South Africa
 Cell: 082 491-7275
 -


   [[alternative HTML version deleted]]

 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo



-- 
Klaus Schliep
Université Paris 6 (Pierre et Marie Curie)
9, Quai Saint-Bernard, 75005 Paris

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] Branch and Bound Maximum Parsimony (Matthew Vavrek)

[Klaus Schliep]

Thanks Liam,

I have two small additions, the first on speeds things up but costs
memory, the second one just simplifies things a bit (and is maybe also
a bit faster).

Regards,
Klaus

On 3/28/11, Liam J. Revell liam.rev...@umb.edu wrote:
 Matthew.

 The only thing that I would add is that if you *really* want to do an
 exhaustive search in R, and your species number is small enough to
 permit an exhaustive search (i.e., =10), then it is straightforward
 enough to do so:

   require(phangorn)
   data-read.phyDat(file=filename,type=USER) # for binary data
   all.trees-allTrees(n=length(data),tip.label=names(data),rooted=FALSE)
all.trees = .uncompressTipLabel(all.trees)
# I know this is ugly
   pscores-vector()
   for(i in 1:length(all.trees))
 pscores[i]-parsimony(all.trees[[i]],data)
# forget the last 3 lines just use
pscores -parsimony(all.trees,data)
# pscores allows multiPhylo objects
   minscore-min(pscores); mp.tree-all.trees[pscores==minscore]

 mp.tree will be a single MP tree or a list if there are several MP trees.

 Of course this will take a long time for more than a 7 or 8 species, and
 will not work at all for more than 10 species.  It's also possible that
 an exhaustive search in a traditional phylogeny inference program like
 PAUP* might be faster if, for instance, PAUP* retains the score for
 parts of the tree that don't change among a set of trees - instead of
 recalculating it each time anew.  That I don't know.

 - Liam

 --
 Liam J. Revell
 University of Massachusetts Boston
 web: http://faculty.umb.edu/liam.revell/
 email: liam.rev...@umb.edu
 blog: http://phytools.blogspot.com

 On 3/28/2011 10:58 AM, Ross Mounce wrote:
 Dear Matthew,


 Branch and Bound searching is often unneccessary, are you sure you need to
 do this? With enough random addition sequences (relative to dataset size)
 you *will* find all MPTs.

 PAUP* despite being familiar to many of us and very fully-featured, is
 definitely ungainly for modern analyses with it's lack of parallelisation.

 Have you considered using TNT? http://www.cladistics.com/aboutTNT.html
 Wiki Manual here: http://tnt.insectmuseum.org/index.php/Main_Page

 Many paleontologists are starting to use this - computationally it's far
 more efficient.
 Just remember, as with any and all phylogenetic analyses it's garbage in,
 garbage out: you must know what you're doing and why *before* you start
 your analysis.


 Kind Regards and good luck,


 Ross Mounce



 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo



-- 
Klaus Schliep
Université Paris 6 (Pierre et Marie Curie)
9, Quai Saint-Bernard, 75005 Paris

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] sort trees by constraint?

[Klaus Schliep]

Dear Brad,

here are some solutions for your problems. I assume unrooted trees
here and ignore the position of the root, so you may additionaly check
for the position of the root.


On 3/24/11, Ruhfel Brad ruh...@gmail.com wrote:
 Dear all,

 I am wondering if there is some way using R to sort a file of multiple trees
 (i.e., BEAST output) to either

 1) return all trees that match a constraint topology (backbone or full
 topology) or

Assume you have all the trees in given in a list of trees (phylo
object) or as multiPhylo object (ape package) and tree0 is the tree
you want them to be compared with:

library(ape)
library(phangorn)
# generate some toy trees
trees = rmtree(10, 5, rooted = FALSE)
tree0 = trees[[5]]
#

tmp = sapply(trees, RF.dist, tree0)  # compute Robinson Foulds distance
result1 = trees[tmp==0]  # the trees you want


 2) return all trees that have a monophyletic clade of a list of taxa.

# Define your partition, here (t1,t2) go together
dat = matrix(c(1,1,0,0,0),dimnames=list(c(t1, t2, t3, t4,
t5),NULL) , ncol=1)
dat = phyDat(dat, type=USER, levels = c(0,1), ambiguity=NULL)

tmp2 = parsimony(trees, dat)
result2 = trees[tmp2==1]  #



 Any help much appreciated!

 Just joined the list... looking forward to it!

 -Brad
 PhD Student
 Harvard Herbaria

 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo



Regards,
Klaus


-- 
Klaus Schliep
Université Paris 6 (Pierre et Marie Curie)
9, Quai Saint-Bernard, 75005 Paris

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] bootstrapping with boot.phylo()

[Klaus Schliep]

Dear Gontran,

for short sequences such a result can happen. Maybe you check whether
any values in
dist.dna(seq100b, model = raw)
are close to 0.75.
I suspect that some elements in this distance matrix are close to 0.75
than pairwise distances for the bootstrap samples are likely to be
equal or greater than 0.75. Most models (K80, JC69 etc.) are not
defined for distances =0.75 and will return Inf or NaN (the 0.75 can
vary a bit, depending on the substitution model). bionj of course does
not like building trees from infinite values as input. With short
sequences the variances are of course larger and you are more likely
to observe this, that's why your larger data set works fine.
However in this cases NaN or Inf are the correct results!
It would be nice to catch this error with a try and use only trees
from finite distance matrices or set infinite values to a large value.
But one should return a warning as these samples are likely to be
biased.

Kind regards,
Klaus




On 3/15/11, Gontran Sonet gontran.so...@naturalsciences.be wrote:

 Dear all,
 When I make a bootstrapping analysis of relatively short DNA sequences
 (100bp), I have an error message:
 Error in bionj(dist.dna(seq100b, model = K80, pairwise.deletion =
 FALSE,  :
NA/NaN/Inf in foreign function call (arg 1)

 I am not able to solve the problem since:
 1)This short sequence data set is part of a larger data set where
 bootstrapping works perfectly using the same commands:

 mytre-bionj(dist.dna(seq100b,model=K80,pairwise.deletion=FALSE,as.matrix=TRUE))
  boot.phylo(mytre,seq100b,FUN=function(seq100b) {
 bionj(dist.dna(seq100b,model=K80,pairwise.deletion=FALSE,
 as.matrix=TRUE))})

 2)The problematic data set doesn't contain any - or ? or N, or
 ambiguities and the dist.dna() doesn't produce any NA or missing value.
seq100b
 396 DNA sequences in binary format stored in a matrix (.
 All sequences of same length: 100
 Labels: XX, XX, XX
 Base composition:
  a c g t
 0.224 0.281 0.222 0.273

 3)This data set also works in other programs (like Mega) without any
 problem.

 Thank you very much for your help

 Gontran Sonet

 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo



-- 
Klaus Schliep
Université Paris 6 (Pierre et Marie Curie)
9, Quai Saint-Bernard, 75005 Paris

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] bind.tree() error

[Klaus Schliep]

Hi David,

have you checked the ordering of the trees? allTrees returns trees in
''pruningwise ordering.   Different functions assume different things
how this should look like, but often don't say it explicitely. Things
like reorder from ape or phangorn:::reorderPruning (a bit faster) are
generally a much faster alternative for read.tree(write.tree(...)),
even this got also faster recently.

Cheers,
Klaus

On 3/8/11, David Bapst dwba...@uchicago.edu wrote:
 Hello all,

 Apologies for sending another post this week, but I've run into a
 strange error with bind.tree(), as below. I'm binding one tree
 (clade) to a tip (4, in this case) on another tree (chosen_tree).
 chosen_tree is a tree from an allTrees() output.

 bind.tree(chosen_tree,clade,where=tip)
 Error in x$edge[sndcol, 2] - newNb[-x$edge[sndcol, 2]] - (n + 2):(n +  :
   number of items to replace is not a multiple of replacement length
 In addition: Warning message:
 In newNb[-x$edge[sndcol, 2]] - (n + 2):(n + x$Nnode) :
   number of items to replace is not a multiple of replacement length

 Bizarrely, while trying to create sample files for you to test, I
 found that write.tree() and read.tree() actually get rid of the
 problem, and furthermore, that the problem is with chosen_tree and
 not clade.  So, as of the moment, I'm unfortunately unable to make
 the error reproducible. Not all trees from an allTrees output do this,
 either.

 bind.tree(chosen_tree,read.tree(text=write.tree(clade,file=)),where=tip)
 Error in x$edge[sndcol, 2] - newNb[-x$edge[sndcol, 2]] - (n + 2):(n +  :
   number of items to replace is not a multiple of replacement length
 In addition: Warning message:
 In newNb[-x$edge[sndcol, 2]] - (n + 2):(n + x$Nnode) :
   number of items to replace is not a multiple of replacement length
 bind.tree(read.tree(text=write.tree(chosen_tree,file=)),clade,where=tip)
 Phylogenetic tree with 5 tips and 4 internal nodes.
 Tip labels:
 [1] 1781   23 24 GR1210 GR3799
 Rooted; no branch lengths.

 Anyone have any idea what's going on? For now, I'll solve it by using
 the read.tree(write.tree()) kludge.
 -Dave


 --
 David Bapst
 Dept of Geophysical Sciences
 University of Chicago
 5734 S. Ellis
 Chicago, IL 60637
 http://home.uchicago.edu/~dwbapst/

 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo



-- 
Klaus Schliep
Université Paris 6 (Pierre et Marie Curie)
9, Quai Saint-Bernard, 75005 Paris

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] The Speed of Different Methods For Obtaining the Descendant Tip Taxa per Node

[Klaus Schliep]

Hi Rob,

this is partly also possible in R. You can find all generic functions
specificly written for an (S3) object using:
methods(class=phylo)
This will miss functions, which are not generic and generic function
which use the default mechanism, e.g. class pml has a specific
logLik.pml  method, but one can use the AIC.default method on
objects of this class.

Cheers,
Klaus




On 3/7/11, Rob Lanfear rob.lanf...@gmail.com wrote:
 Hi All,

 On the issue of visibility of functions in R, I definitely find it a
 struggle to 'discover' things that I could do with a particular object.
 Python has this worked out in a really nice way - if I have an object I can
 simply put a full-stop after it and hit tab, and I get a list of all
 currently loaded functions which could apply to that object. Does R have
 anything similar to this? For instance, if I have a phylo object, is there
 any straightforward way to see which functions currently loaded in R can
 take that object as input? If not, would it be hard to write a function to
 do this?

 Often I find myself reading through the full help pdf for a given of
 libraries to see if I can find what I'm looking for. Depending on what
 people have named their functions this can sometimes be a painful process
 too!

 Rob

 On 8 March 2011 03:35, David Bapst dwba...@uchicago.edu wrote:

 Emmanuel, all-
 That's pretty fast, even on my cheap laptop!

  system.time(desc1-c(as.list(1:Ntip(res_tree)),prop.part(res_tree)))
   user  system elapsed
   0.650.000.66

 As far as why I didn't try prop.part, to be honest, I had no idea that
 prop.part() did that. The help file says:

 Description
 These functions analyse bipartitions found in a series of trees.
 prop.part counts the number of bipartitions found in a series of trees
 given as 

 and

 Value
 prop.part returns an object of class prop.part which is a list with
 an attribute number. The elements of this list are the observed
 clades, and the attribute their respective numbers. If the default
 check.labels = FALSE is used, an attribute labels is added, and the
 vectors of the returned object contains the indices of these labels
 instead of the labels themselves.

 That doesn't suggest to me that it gives the tip descendants of every
 internal node, although that is part of what it does. Until you
 suggested I run it, I had no idea it could give the information I
 wanted. Perhaps a little bit of clarifying of the help file,
 particularly the specific bits of the output, would go a long way.

 Cheers!
 -Dave



 On Sun, Mar 6, 2011 at 10:11 PM, Emmanuel Paradis
 emmanuel.para...@ird.fr wrote:
  Hi David and others,
 
  prop.part() with a single tree does what you want:
 
  set.seed(123)
  res_tree - rtree(1700)
  system.time(desc2 - prop.part(res_tree))
 
user  system elapsed
   0.050   0.000   0.053
 
  edge_end - unique(res_tree$edge[,1])
  system.time(desc1 - Descendants(res_tree,edge_end))
 
user  system elapsed
   0.400   0.000   0.406
 
  This gives the same answer than Descendants:
 
  for (i in seq_along(desc1))
 
  + stopifnot(all(sort(desc1[[i]]) == sort(desc2[[i]])))
 
 
  Klaus makes an excellent reminder about R programming: use vectorisation
  whenever possible (you'll almost always win to do it).
 
  Your examples are very instructive. One lesson is that it is quite easy
  to write R code that tends to be slow. Another one is that it is also
  easy to write different functions doing the same thing. I'm quite
  impressed to see that several packages have re-implemented something
  that has been in ape for some time.
 
  I wish developpers (and users as well) could use these tools from ape
  which have been developed specifically for this kind of objectives
  (fast,
  efficient, and reliable). Is it an issue of 'visibility' of these
 functions?
  This could be a project idea for the next GSoC.
 
  Cheers,
 
  Emmanuel
 
  PS: if you want to speed-up your computations and you have a multi-core
  processor on your machine (which are common on laptops nowadays), you
 could
  try the multicore package (phangorn already supports it):
 
  library(multicore)
  job - parallel(prop.part(res_tree))
  system.time(out - collect(job))
user  system elapsed
   0.080   0.010   0.002
  out - out[[1]]
  identical(out, desc2)
  [1] TRUE
 
  This might be even more interesting with lists of trees.
 
 
  David Bapst wrote on 07/03/2011 04:32:
 
  Klaus-
  Oh, that worked rather splendid! Thanks for letting me know.
 
  system.time(desc-Descendants(res_tree,edge_end))
 
user  system elapsed
1.560.001.56
 
  -Dave
 
  On Sun, Mar 6, 2011 at 3:22 PM, Klaus Schliep klaus.schl...@gmail.com
  wrote:
 
  Hi David,
 
  you can supply Descendents (from phangorn) with a vector instead of
  using sapply. I should have mentioned this in the help file.
  If your vector (edge_end) is long enough (I don't have a exact number
  here) this can be much faster than using sapply as some results are
  reused.
  Here

Re: [R-sig-phylo] Error after ladderizing a tree

[Klaus Schliep]

Hi Alastair,

first here is a simple solution to your other problem:

library(ape)
library(phangorn)
i - 1
#set.seed(i)
r.tree - rtree(20)
#set.seed(i)
sim.data - simSeq(r.tree, l = 200, rootseq=rep('a',200), type=DNA,rate=0.5)
sim.tree - NJ(dist.hamming(sim.data))
sim.tree - optim.parsimony(sim.tree, sim.data, method=sankoff, trace=0)
sim.tree.lad - ladderize(reorder(sim.tree))
plot(sim.tree)
plot(sim.tree.lad) # Gives an error  suggests to use collapse.singles
sim.tree.lad2 - collapse.singles(sim.tree.lad) #returns a tree
plot(sim.tree.lad2)# Gives an error

It seems ladderize expects cladewise ordering.

Cheers,
Klaus


On 3/3/11, Alastair Potts pott...@gmail.com wrote:
 Hi All,
 The code below gives me an error after I ladderize a tree. Any ideas as
 to why?

 library(ape)
 library(phangorn)
 i - 1
 #set.seed(i)
 r.tree - rtree(20)
 #set.seed(i)
 sim.data - simSeq(r.tree, l = 200, rootseq=rep('a',200),
 type=DNA,rate=0.5)
 sim.tree - NJ(dist.hamming(sim.data))
 sim.tree - optim.parsimony(sim.tree, sim.data, method=sankoff, trace=0)
 sim.tree.lad - ladderize(sim.tree)
 plot(sim.tree)
 plot(sim.tree.lad) # Gives an error  suggests to use collapse.singles
 sim.tree.lad2 - collapse.singles(sim.tree.lad) #returns a tree
 plot(sim.tree.lad2)# Gives an error

 Thanks for your time.
 Alastair

 p.s. thanks very much to Joseph for his quick response to my earlier
 query regarding a statistic to compare two trees. His answer of the
 treedist() and RF.dist() functions from the phangorn library were
 exactly what I was looking for.

 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo



-- 
Klaus Schliep
Université Paris 6 (Pierre et Marie Curie)
9, Quai Saint-Bernard, 75005 Paris

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] Problems with transferring trees back into ape from Mesquite

[Klaus Schliep]

Hi David,

there are a few superfluous brackets in your tree (161) and (185)
around single tips. If you remove them read.tree works fine. These
also explains un-collapsed double node branches, you produced them
with adding the brackets.

Cheers,
Klaus

PS: There is a function collapse.singles in ape.

On 3/3/11, David Bapst dwba...@uchicago.edu wrote:
 Hello all,

 I recently decided to manually edit a very large supertree (1700
 taxa) in Mesquite (I had to collapse a few nodes into polytomies,
 something I can only seem to do in Mesquite). I read it out of R as a
 Newick file, than into Mesquite, where I did my little edits and then
 saved the tree as a Nexus file, which I then tried to read into ape.
 But then I get this.

 spec_tree-read.nexus(file.choose())
 Error in if (sum(tr$edge[, 1] == ROOT) == 1  dim(tr$edge)[1]  1) { :
   missing value where TRUE/FALSE needed

 I also tried reading the Nexus file into TreeFig, which didn't work.
 If I cut out a lot of the internal structure of the Nexus file (who
 needs Mesquite and Taxon blocks anyway?), I can get TreeFig to read
 it, export it into a Newick file or a new Nexus file and then try
 opening it in R with read.nexus(). That works. But then i can't write
 it in a new file in R with write.tree() nor write.nexus(), and also
 plot(tree) causes R to freeze.

 write.tree(spec_id_tree,file.choose())
 Error in kids[[parent[i]]] : attempt to select less than one element

 write.nexus(spec_id_tree,file.choose())
 #NEXUS
 [R-package APE, Mon Feb 28 10:03:59 2011]

 Read 1756 items
 Error in 1:(start - 1) : argument of length 0

 After saving and opening the file in Mesquite a bunch of times, I
 repeatedly found new un-collapsed double node branches on the tree. No
 idea where they came from (maybe a by-product of my node collapsing?).
 I removed the new ones that popped up, deassigned branch lengths to
 the tree (it's a super-cladogram) and then saved and... finally, R
 read the nexus file and write.tree() worked. write.nexus() didn't work
 (same error as above), but I can read my tree and write it, at least,
 so my issue is solved but I'd still like to know what was causing the
 problems to begin with. I have no idea where the un-collapsed double
 branches came from or why they were causing this issue.

 I am using R v2.12.2 and ape v2.6-3. So that the error is completely
 reproducible, I have placed two Nexus files on my website which
 produce both the read.nexus() error (trash_a.nex) and the
 write.nexus() and write.tree() errors (trash_b). These are small,
 partially collapsed portions of my supertree with the species names
 deleted. They can be found at:

 http://home.uchicago.edu/~dwbapst/trash_a.nex

 http://home.uchicago.edu/~dwbapst/trash_b.nex

 -Dave


 --
 David Bapst
 Dept of Geophysical Sciences
 University of Chicago
 5734 S. Ellis
 Chicago, IL 60637
 http://home.uchicago.edu/~dwbapst/

 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo



-- 
Klaus Schliep
Université Paris 6 (Pierre et Marie Curie)
9, Quai Saint-Bernard, 75005 Paris

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] Retention Index

[Klaus Schliep]

Hi Miguel,

phangorn can use any discrete set of variables, not only DNA data. I
deleted the (dna) in the help to avoid confusion in future versons of
phangorn.
The RI function however needs an object of class phyDat as input. You
probably need to do something similar like this line to convert your
data:
myData = phyDat(myOldData, type = USER, levels = c(1, 2,3))
vignette(phangorn-specials) gives also some details and hints.

You run these lines, where ryData contains binary Data (the states are
r and y):
library(phangorn)
example(parsimony)
RI(treeRatchet, Laurasiatherian)
(ryData - acgt2ry(Laurasiatherian) )
RI(treeRatchet, ryData)
str(as.data.frame(ryData))


Kind regards
Klaus


On 2/18/11, Miguel Verdu miguel.ve...@uv.es wrote:
 Hi all:

 Is possible to calculate the Retention Index of a
 discrete trait in a phylogenetic tree with R?.
 I´ve tried with phangorn but it seems it only works with DNA data.

 Thanks in advance

 Miguel Verdú





 ***
   Centro de Investigaciones sobre Desertificacion -CIDE-
   (CSIC/UV/GV)
   Apartado Oficial
   46470 ALBAL, VALENCIA (SPAIN)
   phone: (+34) 961220540 ext. 110
   fax: (+34) 961270967
   miguel.ve...@uv.es
 http://www.uv.es/verducam

 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo



-- 
Klaus Schliep
Université Paris 6 (Pierre et Marie Curie)
9, Quai Saint-Bernard, 75005 Paris

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Re: [R-sig-phylo] optimizing branch length of a tree topology based on a distance matrix

[Klaus Schliep]

Hi Sebastian,

the phangorn  package contains a function to construct design matrices
from a tree (class phylo from ape).

library(phangorn)
X = designTree(tree, method = unrooted)
dm = as.matrix(dm)[tree$tip, tree$tip]
y = dm[lower.tri(dm)]
lm(y~X-1) # a normal Least Squares fit

There are several packages in R to solve quadratic programming
problems (non-negative LS in your case), e.g. function solve.QP in the
package quadprog should do the trick. So it should be straight forward
to solve your problem, once you figured out the parameters of
solve.QP.


Regards
Klaus

On 2/17/11, Sébastien Lavergne sebastien.laver...@ujf-grenoble.fr wrote:
 Dear all,

 I would like to optimise branch length of a known tree topology based on
 a distance matrix (non genetic data), and prohibiting negative branch
 lengths.
 Is there a way to do this in R ?
 Any help or hints external to R are also appreciated (I just know there
 is way through this in PAUP,, by constraining 'hsearch' function with a
 tree topology).

 Thanks in advance
 Seb


 --

 -
  Sébastien Lavergne
  Laboratoire d'Ecologie Alpine, UMR-CNRS 5553
  Université Joseph Fourier
  BP 53, 38041 Grenoble Cedex 9, France
  tel +33 (0)4 76 63 54 50
  http://www-leca.ujf-grenoble.fr/membres/lavergne.htm

 ___
 R-sig-phylo mailing list
 R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo



-- 
Klaus Schliep
Université Paris 6 (Pierre et Marie Curie)
9, Quai Saint-Bernard, 75005 Paris

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo