Tuesday, August 23, 2016

Empirical co-phylogenetic ("cophylo") plotting

In some recent updates I have added some functionality to the methods associated with the phytools object class "cophylo", for co-phylogenetic plotting.

I'm considering writing up a short paper describing the methods, and some of the optimization algorithms that I have designed to maximize tip alignment by different criteria, so I thought it might be useful to find an empirical 'exemplar' dataset.

Today I found the following article by Lopez-Vaamonde et al. (2001) from MPE:

Lopez-Vaamonde, C., J. Y. Rasplus, G. D. Weiblen, & J. M. Cook. 2001. Molecular phylogenies of fig wasps: Partial cocladogenesis of pollinators and parasites. Molecular Phylogenetics and Evolution 21: 55-71.

Rather than try to find the trees online, I just transcribed them into Newick cladograms by eye. Here is the result:

library(phytools)
Pleistodontes<-read.tree(text="(((1,(2,3)100)50,(((((((4,(5,(6,7)100)69)100,(8,(9,10))81)61,11),(12,13)100),(14,15)100),(16,17)96),18)),19);")
Pleistodontes$tip.label<-c("P._froggatti",
    "P._blandus",
    "P._near_blandus",
    "P._proximus",
    "P._athysanus",
    "P._cuneatus",
    "P._astrabocheilus",
    "P._macrocainus",
    "P._rigisamos",
    "P._imperialis",
    "P._nigriventris",
    "P._xanthocephalus",
    "P._greenwoodi",
    "P._rieki",
    "P._plebejus",
    "P._schizodontes",
    "P._nitens",
    "P._spec._nov._1",
    "P._regalis")
## rotate all nodes to match orientation in article
Pleistodontes<-rotateNodes(Pleistodontes,"all")
Sycoscapter<-read.tree(text="(1,((((2,3)100,((4,5),(6,7)78)100),((8,9)50,10)99),(((11,12),13)91,(14,15)77)));")
Sycoscapter$tip.label<-c("S._australis_{macrophylla}",
    "S._3_{playpoda}",
    "S._5_{cerasicarpa}",
    "S._6_{lilliputiana}",
    "S._7_{subpuberula}",
    "S._1_{brachypoda}",
    "S._15_{rubiginosa}",
    "S._9_{triradiata}",
    "S._10_{crassipes}",
    "S._8_{pleurocarpa}",
    "S._11_{glandifera}",
    "S._12_{xylosycia}",
    "S._14_{hesperidiiformis}",
    "S._2_{aff._obliqua}",
    "S._4_{obliqua}")
Sycoscapter<-rotateNodes(Sycoscapter,"all")
assoc<-cbind(
    Pleistodontes$tip.label[c(7,8,5,6,18,1,3,4,10,15,13,16,12,14,19)],
    Sycoscapter$tip.label)
assoc

##       [,1]                [,2]                        
##  [1,] "P._greenwoodi"     "S._4_{obliqua}"            
##  [2,] "P._xanthocephalus" "S._2_{aff._obliqua}"       
##  [3,] "P._plebejus"       "S._14_{hesperidiiformis}"  
##  [4,] "P._rieki"          "S._12_{xylosycia}"         
##  [5,] "P._blandus"        "S._11_{glandifera}"        
##  [6,] "P._regalis"        "S._8_{pleurocarpa}"        
##  [7,] "P._nitens"         "S._10_{crassipes}"         
##  [8,] "P._schizodontes"   "S._9_{triradiata}"         
##  [9,] "P._imperialis"     "S._15_{rubiginosa}"        
## [10,] "P._athysanus"      "S._1_{brachypoda}"         
## [11,] "P._astrabocheilus" "S._7_{subpuberula}"        
## [12,] "P._proximus"       "S._6_{lilliputiana}"       
## [13,] "P._macrocainus"    "S._5_{cerasicarpa}"        
## [14,] "P._cuneatus"       "S._3_{playpoda}"           
## [15,] "P._froggatti"      "S._australis_{macrophylla}"

Now let's create our "cophylo" object without rotation:

obj<-cophylo(Pleistodontes,Sycoscapter,assoc=assoc,
    rotate=FALSE)
obj
plot(obj,link.type="curved",link.lwd=3,link.lty="solid",
    link.col=make.transparent("blue",0.25),fsize=0.8)

This matches the published co-phylogenetic plot as follows:

Of course, we can see if we can 'improve' on the tip matching by any criterion.

cophylo allows us to optimize the rotation based on basically any arbitrary criterion.

The default is to minimize the sum of squares difference between the vertical positions of the tips in both trees; however we can also minimize the sum of the absolute values. For instance:

obj<-cophylo(Pleistodontes,Sycoscapter,assoc=assoc,
    print=TRUE)

## Rotating nodes to optimize matching...

## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778
## objective: 359.577777777778

## objective: 357.777777777778
## objective: 357.777777777778

## objective: 345.111111111111
## objective: 345.111111111111
## objective: 345.111111111111

## objective: 322.777777777778

## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856

## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778
## objective: 322.777777777778

## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856
## objective: 201.17728531856

## Done.

plot(obj,link.type="curved",link.lwd=3,link.lty="solid",
    link.col=make.transparent("blue",0.25),fsize=0.8)
## absolute values
obj<-cophylo(Pleistodontes,Sycoscapter,assoc=assoc,
    fn=function(x) abs(x),print=TRUE)

## Rotating nodes to optimize matching...

## objective: 57.2666666666667

## objective: 57.2666666666667
## objective: 57.2666666666667
## objective: 57.2666666666667
## objective: 57.2666666666667
## objective: 57.2666666666667
## objective: 57.2666666666667
## objective: 57.2666666666667
## objective: 57.2666666666667
## objective: 57.2666666666667
## objective: 57.2666666666667
## objective: 57.2666666666667

## objective: 56.2666666666667
## objective: 56.2666666666667

## objective: 52.4
## objective: 52.4
## objective: 52.4

## objective: 51.4

## objective: 40.5789473684211
## objective: 40.5789473684211
## objective: 40.5789473684211
## objective: 40.5789473684211
## objective: 40.5789473684211
## objective: 40.5789473684211
## objective: 40.5789473684211
## objective: 40.5789473684211
## objective: 40.5789473684211
## objective: 40.5789473684211
## objective: 40.5789473684211
## objective: 40.5789473684211
## objective: 40.5789473684211
## objective: 40.5789473684211

## objective: 51.4
## objective: 51.4
## objective: 51.4
## objective: 51.4
## objective: 51.4
## objective: 51.4
## objective: 51.4
## objective: 51.4
## objective: 51.4
## objective: 51.4
## objective: 51.4
## objective: 51.4
## objective: 51.4
## objective: 51.4
## objective: 51.4
## objective: 51.4
## objective: 51.4
## objective: 51.4

## objective: 40.5789473684211
## objective: 40.5789473684211
## objective: 40.5789473684211
## objective: 40.5789473684211
## objective: 40.5789473684211
## objective: 40.5789473684211
## objective: 40.5789473684211
## objective: 40.5789473684211
## objective: 40.5789473684211
## objective: 40.5789473684211
## objective: 40.5789473684211
## objective: 40.5789473684211
## objective: 40.5789473684211
## objective: 40.5789473684211

## Done.

plot(obj,link.type="curved",link.lwd=3,link.lty="solid",
    link.col=make.transparent("blue",0.25),fsize=0.8)

plot of chunk unnamed-chunk-2

It is even possible to, for instance, choose a function that does not make sense to optimize matching - for instance the inverse of the sum of squares:

obj<-cophylo(Pleistodontes,Sycoscapter,assoc=assoc,
    fn=function(x) 1/x^2,print=TRUE)

## Rotating nodes to optimize matching...

## objective: 277.233847021676

## objective: 4.77348935389099
## objective: 4.77348935389099

## objective: 4.75539116602117
## objective: 4.75539116602117
## objective: 4.75539116602117
## objective: 4.75539116602117

## objective: 3.7739829235819
## objective: 3.7739829235819
## objective: 3.7739829235819
## objective: 3.7739829235819
## objective: 3.7739829235819
## objective: 3.7739829235819
## objective: 3.7739829235819
## objective: 3.7739829235819
## objective: 3.7739829235819
## objective: 3.7739829235819
## objective: 3.7739829235819

## objective: 6.05514593516918

## objective: 2.75709265189575

## objective: 1.56956478408719

## objective: 1.55965836490761
## objective: 1.55965836490761
## objective: 1.55965836490761
## objective: 1.55965836490761
## objective: 1.55965836490761
## objective: 1.55965836490761
## objective: 1.55965836490761

## objective: 0.506025805166682
## objective: 0.506025805166682
## objective: 0.506025805166682
## objective: 0.506025805166682

## objective: 0.307384379625682
## objective: 0.307384379625682

## objective: 0.307379556534829
## objective: 0.307379556534829
## objective: 0.307379556534829
## objective: 0.307379556534829
## objective: 0.307379556534829
## objective: 0.307379556534829
## objective: 0.307379556534829
## objective: 0.307379556534829
## objective: 0.307379556534829
## objective: 0.307379556534829

## objective: 0.306253296306206

## objective: 0.303195941167207
## objective: 0.303195941167207

## objective: 0.302741287258059

## objective: 0.301049328645266
## objective: 0.301049328645266

## objective: 0.483016922848627
## objective: 0.483016922848627
## objective: 0.483016922848627
## objective: 0.483016922848627
## objective: 0.483016922848627
## objective: 0.483016922848627
## objective: 0.483016922848627
## objective: 0.483016922848627
## objective: 0.483016922848627
## objective: 0.483016922848627
## objective: 0.483016922848627
## objective: 0.483016922848627
## objective: 0.483016922848627
## objective: 0.483016922848627

## objective: 0.301049328645266
## objective: 0.301049328645266
## objective: 0.301049328645266
## objective: 0.301049328645266
## objective: 0.301049328645266
## objective: 0.301049328645266
## objective: 0.301049328645266
## objective: 0.301049328645266
## objective: 0.301049328645266
## objective: 0.301049328645266
## objective: 0.301049328645266
## objective: 0.301049328645266

## objective: 0.30081492049661
## objective: 0.30081492049661
## objective: 0.30081492049661

## objective: 0.300433102030201
## objective: 0.300433102030201
## objective: 0.300433102030201

## objective: 0.482028221479567
## objective: 0.482028221479567
## objective: 0.482028221479567
## objective: 0.482028221479567
## objective: 0.482028221479567
## objective: 0.482028221479567
## objective: 0.482028221479567
## objective: 0.482028221479567
## objective: 0.482028221479567
## objective: 0.482028221479567
## objective: 0.482028221479567
## objective: 0.482028221479567
## objective: 0.482028221479567
## objective: 0.482028221479567

## objective: 0.300433102030201
## objective: 0.300433102030201
## objective: 0.300433102030201
## objective: 0.300433102030201
## objective: 0.300433102030201
## objective: 0.300433102030201
## objective: 0.300433102030201
## objective: 0.300433102030201
## objective: 0.300433102030201
## objective: 0.300433102030201
## objective: 0.300433102030201
## objective: 0.300433102030201
## objective: 0.300433102030201
## objective: 0.300433102030201
## objective: 0.300433102030201
## objective: 0.300433102030201
## objective: 0.300433102030201
## objective: 0.300433102030201

## objective: 0.482028221479567
## objective: 0.482028221479567
## objective: 0.482028221479567
## objective: 0.482028221479567
## objective: 0.482028221479567
## objective: 0.482028221479567
## objective: 0.482028221479567
## objective: 0.482028221479567
## objective: 0.482028221479567
## objective: 0.482028221479567
## objective: 0.482028221479567
## objective: 0.482028221479567
## objective: 0.482028221479567
## objective: 0.482028221479567

## Done.

plot(obj,link.type="curved",link.lwd=3,link.lty="solid",
    link.col=make.transparent("blue",0.25),fsize=0.8)

plot of chunk unnamed-chunk-3

This will have the tendency to minimize matching - which, of course, we don't want.

Finally, we can obviously add annotation - such as the bootstrap values at nodes from the original paper. (Though to get the look we want can be tricky.) For instance:

obj<-cophylo(Pleistodontes,Sycoscapter,assoc=assoc,
    rotate=FALSE)
plot(obj,link.type="curved",link.lwd=3,link.lty="solid",
    link.col=make.transparent("blue",0.25),fsize=0.8)
nodes<-1:obj$trees[[1]]$Nnode+Ntip(obj$trees[[1]])
nodelabels.cophylo(node=c(30,33),
    text=obj$trees[[1]]$node.label[c(30,33)-Ntip(obj$trees[[1]])],
    adj=c(1.1,-0.2),
    frame="none",cex=0.7,which="left")
nodelabels.cophylo(node=35,
    text=obj$trees[[1]]$node.label[35-Ntip(obj$trees[[1]])],
    adj=c(-0.1,-0.1),
    frame="none",cex=0.7,which="left")
nodelabels.cophylo(node=setdiff(nodes,c(30,33,35)),
    text=obj$trees[[1]]$node.label[
    setdiff(nodes,c(30,33,35))-Ntip(obj$trees[[1]])],
    adj=c(1.1,1.1),
    frame="none",cex=0.7,which="left")
nodes<-1:obj$trees[[2]]$Nnode+Ntip(obj$trees[[2]])
nodelabels.cophylo(node=c(19,20,24),
    text=obj$trees[[2]]$node.label[c(19,20,24)-Ntip(obj$trees[[2]])],
    adj=c(-0.2,-0.4),
    frame="none",cex=0.7,which="right")
nodelabels.cophylo(node=setdiff(nodes,c(19,20,24)),
    text=obj$trees[[2]]$node.label[setdiff(nodes,c(19,20,24))-
    Ntip(obj$trees[[2]])],
    adj=c(-0.2,1.1),
    frame="none",cex=0.7,which="right")

plot of chunk unnamed-chunk-4

That's it.

Posted by at

5 comments:

  1. dwbapstAugust 24, 2016 at 10:10 AM

    I think having a citation for us to use for your cophylo method would be good; at the moment I've been just citing phytools in some submitted manuscripts.

    Have you given any time toward better optimization for rotating trees with polytomies? I've considered doing this myself, since I always seem to be applying cophylo to trees with big polytomies...

    ReplyDelete
    Replies
    1. Liam RevellAugust 24, 2016 at 1:56 PM

      Dave. It's funny that you say that. I just started working on a short ms. - probably to submit to MEE.

      With regards to your second question - yes in fact! See here: http://blog.phytools.org/2016/08/co-phylogenetic-plotting-for-trees-with.html.

      Thanks for the comments Dave. - Liam

      Delete
    2. dwbapstAugust 25, 2016 at 6:11 PM

      1) Liam, yeah, you mentioned that you were working on such a manuscript in this post. I was agreeing that it was a good idea.

      2) Ah, I missed this while I was off where I had no internet for a week. Cool! I'll have to try the multifurcation rotating out... when some manuscripts return from review.

      Cheers,
      -Dave

      Delete
  2. David KnMarch 10, 2022 at 8:59 AM

    Dear Liam,

    Is there a way to display boostrap values to a cophylogenetic tree, in panel A and B, capped at a certain threshold (say 80% bootstrap), and doing this without invoking separate "nodelabels.cophylo(node=c(19,20,24), ..." commands?

    Thanks in advance for the support.
    Bw,
    David

    ReplyDelete
  3. ElenaSeptember 7, 2023 at 8:12 AM

    Thank you very much Liam, really useful!
    Quick question, is there a way to save the cophylo in the environment as a whole? (with the node labels and coloured lines etc) or is there a way to transform it to phylo so we can add more data to the cophylo file?
    Thank you very much in advance!!

    ReplyDelete

Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.

Subscribe to: Post Comments (Atom)