Monday, November 2, 2015

phylo.to.map with right-facing trees; preliminary t-shirt design for Valdivia workshop

I will be teaching a workshop at Universidad Austral de Chile in Valdivia in December. Right now I'm working on the design for the t-shirt that we have made a tradition in prior similar workshops (e.g., 1, 2).

I decided it might be cool to project a phylogeny onto the geographic outline of Chile. The only problem is that Chile is a very long country, and phylo.to.map previously only permits a projection of the tree from above facing down.

Well - that has now been addressed. I just pushed a phytools update that allows a phylogeny to be plotted to the left of a geographic map facing rightwards. Here is a demo.

First, we need to install from GitHub:

library(devtools)
install_github("liamrevell/phytools")

Load packages:

library(phytools)
library(phangorn) ## we will use phangorn too

For fun, let's use the lat. & long. coordinates of some major Chilean cities:

X<-read.csv("Chile-cities.csv",row.names=1)
X<-as.matrix(X)
X
##                 Lat   Long
## Santiago     -33.45 -70.67
## Valparaiso   -33.05 -71.62
## Antofagasta  -23.65 -70.40
## Temuco       -38.75 -72.67
## Concepcion   -36.82 -73.05
## Rancagua     -34.17 -70.75
## Talca        -35.43 -71.67
## Arica        -18.48 -70.33
## Iquique      -20.22 -70.15
## Puerto_Montt -41.47 -72.93
## La_Serena    -29.90 -71.25
## Chillan      -36.60 -72.12
## Osorno       -40.57 -73.14
## Valdivia     -39.75 -72.50
## Calama       -22.47 -68.93
## Copiapo      -27.36 -70.33
## Los_Angeles  -37.47 -72.35
## Punta_Arenas -53.17 -70.93

(Hopefully I got most of these right!)

Rather than use a random “tree” - I figured why not compute a tree from the latitudinal distances between cities:

D<-dist(X[,1])
D
##              Santiago Valparaiso Antofagasta Temuco Concepcion Rancagua
## Valparaiso       0.40                                                  
## Antofagasta      9.80       9.40                                       
## Temuco           5.30       5.70       15.10                           
## Concepcion       3.37       3.77       13.17   1.93                    
## Rancagua         0.72       1.12       10.52   4.58       2.65         
## Talca            1.98       2.38       11.78   3.32       1.39     1.26
## Arica           14.97      14.57        5.17  20.27      18.34    15.69
## Iquique         13.23      12.83        3.43  18.53      16.60    13.95
## Puerto_Montt     8.02       8.42       17.82   2.72       4.65     7.30
## La_Serena        3.55       3.15        6.25   8.85       6.92     4.27
## Chillan          3.15       3.55       12.95   2.15       0.22     2.43
## Osorno           7.12       7.52       16.92   1.82       3.75     6.40
## Valdivia         6.30       6.70       16.10   1.00       2.93     5.58
## Calama          10.98      10.58        1.18  16.28      14.35    11.70
## Copiapo          6.09       5.69        3.71  11.39       9.46     6.81
## Los_Angeles      4.02       4.42       13.82   1.28       0.65     3.30
## Punta_Arenas    19.72      20.12       29.52  14.42      16.35    19.00
##              Talca Arica Iquique Puerto_Montt La_Serena Chillan Osorno
## Valparaiso                                                            
## Antofagasta                                                           
## Temuco                                                                
## Concepcion                                                            
## Rancagua                                                              
## Talca                                                                 
## Arica        16.95                                                    
## Iquique      15.21  1.74                                              
## Puerto_Montt  6.04 22.99   21.25                                      
## La_Serena     5.53 11.42    9.68        11.57                         
## Chillan       1.17 18.12   16.38         4.87      6.70               
## Osorno        5.14 22.09   20.35         0.90     10.67    3.97       
## Valdivia      4.32 21.27   19.53         1.72      9.85    3.15   0.82
## Calama       12.96  3.99    2.25        19.00      7.43   14.13  18.10
## Copiapo       8.07  8.88    7.14        14.11      2.54    9.24  13.21
## Los_Angeles   2.04 18.99   17.25         4.00      7.57    0.87   3.10
## Punta_Arenas 17.74 34.69   32.95        11.70     23.27   16.57  12.60
##              Valdivia Calama Copiapo Los_Angeles
## Valparaiso                                      
## Antofagasta                                     
## Temuco                                          
## Concepcion                                      
## Rancagua                                        
## Talca                                           
## Arica                                           
## Iquique                                         
## Puerto_Montt                                    
## La_Serena                                       
## Chillan                                         
## Osorno                                          
## Valdivia                                        
## Calama          17.28                           
## Copiapo         12.39   4.89                    
## Los_Angeles      2.28  15.00   10.11            
## Punta_Arenas    13.42  30.70   25.81       15.70
tree<-untangle(upgma(D),"read.tree") ## untangle is hack

(Use of untangle is a hack to ensure that the tree is in good conformation for phylo.to.map.)

Now let's proceed!

obj<-phylo.to.map(tree,X,regions="Chile",direction="rightwards",
    plot=FALSE)
## objective: 32
## objective: 32
## objective: 32
## objective: 26
## objective: 26
## objective: 26
## objective: 8
## objective: 2
## objective: 2
## objective: 2
## objective: 2
## objective: 2
## objective: 2
## objective: 2
## objective: 2
## objective: 2
## objective: 0
obj
## Object of class "phylo.to.map" containing:
## 
## (1) A phylogenetic tree with 18 tips and 17 internal nodes.
## 
## (2) A geographic map with range:
##      -55.89N, -17.51N
##      -109.43W, -66.44W.
## 
## (3) A table containing 18 geographic coordinates.
plot(obj,direction="rightwards",ftype="i")

plot of chunk unnamed-chunk-5

One issue here is the spacing - which actually turns out to be the result of the inclusion in “Chile” of Chile's pacific island territories, such as Easter Island. Let's unceremoniously chop off those islands!

ii<-which(obj$map$x<(-77))
obj$map$x<-obj$map$x[-ii]
obj$map$y<-obj$map$y[-ii]
obj$map$range[1]<--81
plot(obj,direction="rightwards",split=c(0.6,0.4),ftype="i")

plot of chunk unnamed-chunk-6

Finally, let's plot it with the tip labels removed & a special star at the location of our upcoming workshop!

plot(obj,direction="rightwards",split=c(0.3,0.7),ftype="off",
    colors=c(grey(0.3),"black"))
points(X["Valdivia",2],X["Valdivia",1],pch=8,cex=1.2)

plot of chunk unnamed-chunk-7

Not bad! Of course, as always, we can export a higher quality version using pdf:

pdf(file="phylo.to.map-Chile.pdf",width=5,height=8)
plot(obj,direction="rightwards",split=c(0.3,0.7),ftype="off",
    colors=c(grey(0.3),"black"))
points(X["Valdivia",2],X["Valdivia",1],pch=8,cex=1.2)
dev.off()
## windows 
##       2

Here it is!

That's all for now.

10 comments:

  1. Excellent spatial data representation!! I'd like to do the same for my data, but instead of to use upgma method for inferring phylogenetic relationship, can I use any Bayesian method? How?
    Thanks in advance
    Regards
    Carlos

    ReplyDelete
    Replies
    1. Hola Carlos. Usé UPGMA solamente para generar esta demostración del método. Se puede utilizar un árbol estimado por algún método que quieres, incluyendo la inferencia bayesiana. Si quieres representar datos espaciales como dendrogramo, en lugar de una filogenia evolutiva, es posible de convertir entre objetos de clase “hclust” y objetos de clase “phylo” (que necesita para usar ésta función) usando la función as.phylo en ape. ¡Déjame saber si tienes algunas dudas o preguntas sobre este tema! - Liam

      Delete
    2. This comment has been removed by the author.

      Delete
    3. Gracias por tu rápida respuesta. Creo entender lo que dices, pero debo probar para ver como funciona. Por otro lado me gustaría saber sobre el curso o workshop que impartirás en Valdivia. Es abierto para extranjeros? Dónde me podría inscribir?
      Estamos en contacto!!
      Carlos

      Delete
    4. Hola Carlos.

      El taller que estamos haciendo en Valdivia en diciembre fue organizado por Dr. Roberto Nespolo en la Universidad Austral de Chile (http://www.nespolo.cl/). Si tienes interés en este curso, tienes que contactar a Roberto, aunque no estoy seguro si este taller en particular es abierto a gente afuera de UACh ni si hay todavía espacio por estudiantes adicionales.

      Estoy también enseñando un curso de duración similar, pero en inglés, en marzo con Transmitting Science. Hay más información sobre este taller en la página de los organizadores: http://www.transmittingscience.org/courses/evol/phytools/.

      Finalmente, estoy actualmente pidiendo fondos para ofrecer más versiones de mi curso en métodos comparados en castellano en varios sitios latinoamericanos. Entonces, aunque todavía tenemos que cruzar los dedos, hay una posibilidad fuerte de varias ediciones de más en los años que vienen.Obviamente, algunos cursos abiertos a inscripción pública serían anunciados aquí en este blog entre otros sitios.

      Gracias por tu interes en phytools y por los comentarios! - Liam

      Delete
  2. Estimado Liam,
    Gracias por tu pronta respuesta. Por la cercanía a mi país (Bolivia), me gustaría participar en el curso taller de Valdivia. Escribiré a Roberto, haber que dice. Caso contrarío tengo "keep cross fingers" y esperar que tus fondos tengan éxito para que puedas venir a Bolivia o algún otro vecino país.
    Suerte y mantenemos el contacto.
    Carlos

    ReplyDelete
    Replies
    1. Buenos días Carlos.
      Si no logras en asistir en Valdivia, debes saber que si estamos exitosos en obtener fondos para apoyar versiones adicionales del taller, tendremos también becas para apoyar viaje y alojamiento por los participantes. Esperamos ofrecer a menos uno o dos de estos talleres en Sudamérica, y uno o dos más en el Caribe o en América Central. La inscripción cada vez sería limitada a veinte-cinco individuos, más o menos.
      Saludos, Liam

      Delete
  3. Liam, thank you very much for this post. It has been very helpful and I already was able to build and plot the phymap object for my data. However, I was wondering if it is possible to add other spatial features to the map. For example, is it possible to add raster layers (e.g. to show elevation)? Any insight and help is very much appreciated. Thank you in advance.
    Best,
    Laura Céspedes

    ReplyDelete
    Replies
    1. I would be interested in this feature as well!

      Delete
  4. Whenever I try to plot I get objective = 0 and the following error

    Error in coords[cw$tip.label, 2:1] : subscript out of bounds

    ReplyDelete

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