Python for Bioinformatics: Progress with RPy

Friday, November 19, 2010

Progress with RPy

I made significant progress this morning on using RPy. I was having trouble calling a function from the ape package in R. The problem was simply a matter of namespaces, or "environments" as R calls them (I know they're more complex than that but ..).

So if I want the function read.tree, I have to specify it's from ape and also the syntax is different (standard call first, ape call second):

from rpy2 import robjects
robjects.r['pi'][0]
from rpy2.robjects.packages import importr
ape = importr('ape')
read_tree = ape.read_tree
tr = read_tree('/Users/telliott_admin/Desktop/tree.txt')

>>> tr

[57467 refs]

>>> print tr

Phylogenetic tree with 6 tips and 4 internal nodes.

Tip labels:
[1] "Stenotrophomonas_maltophilia" "Kingella_oralis"             
[3] "Pseudomonas_aeruginosa"       "Salmonella_typhi"            
[5] "Escherichia_coli"             "Haemophilus_parainfluenzae"  

Unrooted; includes branch lengths.

[57468 refs]

>>> for e in tr:
...     print e
... 
      [,1] [,2]
 [1,]    7    8
 [2,]    8    1
 [3,]    8    2
 [4,]    7    3
 [5,]    7    9
 [6,]    9   10
 [7,]   10    4
 [8,]   10    5
 [9,]    9    6

[1] 4

[1] "Stenotrophomonas_maltophilia" "Kingella_oralis"             
[3] "Pseudomonas_aeruginosa"       "Salmonella_typhi"            
[5] "Escherichia_coli"             "Haemophilus_parainfluenzae"  

[1] 0.00827 0.07574 0.08026 0.05950 0.03863 0.03356 0.01297 0.01491 0.06113

[57729 refs]

>>> list(tr.names)
['edge', 'Nnode', 'tip.label', 'edge.length']
[57626 refs]

The second thing is that I would normally call plot with a tree and an argument tip.color= ...

plot(tr, tip_color="black")
Warning messages:
1: In plot.window(...) : "tip_color" is not a graphical parameter
2: In plot.xy(xy, type, ...) : "tip_color" is not a graphical parameter

[57936 refs]

Instead of doing that, what worked for me is to call the specific plot function from ape:

>>> print ape.plot_phylo
function (x, type = "phylogram", use.edge.length = TRUE, node.pos = NULL, 
    show.tip.label = TRUE, show.node.label = FALSE, edge.color = "black", 
    edge.width = 1, edge.lty = 1, font = 3, cex = par("cex"), 
    adj = NULL, srt = 0, no.margin = FALSE, root.edge = FALSE, 
    label.offset = 0, underscore = FALSE, x.lim = NULL, y.lim = NULL, 
    direction = "rightwards", lab4ut = "horizontal", tip.color = "black", 
    ...) 

ape.plot_phylo(tr, tip_color='black')

Next up, trying not to deal with the plotting window. RPy is set up to use X11, and it doesn't look so hot. So, I'd like to save the plots as pdf files. I found the instructions here:


>>> grdevices.pdf('/Users/telliott_admin/Desktop/plot.pdf')
<RObject - Python:0x26bd9f8 / R:0x8bdfb0>
[59806 refs]
>>> ape.plot_phylo(tr, tip_color='red', cex=2)
<Vector - Python:0x26b2b78 / R:0x2d57c30>
[59825 refs]
>>> grdevices.dev_off()
<IntVector - Python:0x26bdbb8 / R:0xd746a8>
[59845 refs]

The other issue I really haven't solved is RPy crashing. Or rather, I can't reproduce it with Python 2.6.6 installed. I did this by building it with debug enabled (That's where the lines like [57936 refs]are coming from in the output above. But that's another story.