Wednesday, February 26, 2014

Rphylip mostly done....

I've been continuing to plug away at Rphylip, our new R interface for Joe Felsenstein's phylogeny methods package PHYLIP. More information about this effort can be found here: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10. The latest additions are Rkitsch, an interface for KITSCH (a least-squares & ME inference program with a clock constraint); Rrestdist, an interface for RESTDIST (a program for distance calculation from restriction site or fragment data); Rrestml, an interface for RESTML (a restriction site ML tree inference program); and Rclique, an interface for CLIQUE (phylogeny inference from binary characters by the compatibility method). I also added a new object class, "rest.data", for restriction data, as well as some methods to convert to the class & print.

Rphylip is on GitHub, and the latest build can be downloaded & installed from source.

Here's a quick demo of Rrestml using the demo dataset from the RESTML documentation in PHYLIP:

> require(Rphylip)
Loading required package: Rphylip
Loading required package: ape
> packageVersion("Rphylip")
[1] ‘0.1.20’
> data(restriction.data)
> restriction.data
13 restriction site scores for 5 species stored in a object of class "rest.data".

All sequences of same length: 13

Number of restriction enzymes used to generate the data: 2

Labels: Alpha Beta Gamma Delta Epsilon

> mltree<-Rrestml(restriction.data,quiet=TRUE)

Restriction site Maximum Likelihood method, version 3.695


  Recognition sequences all 6 bases long

Sites absent from all species are assumed to have been omitted


     +2        
  +--1 
  |  |  +4        
  |  +--2 
  |     +5        
  | 
  3----3        
  | 
  +1        


remember: this is an unrooted tree!

Ln Likelihood =   -40.31850


Between  And  Length   Approx. Confidence Limits
-------  ---  ------   ------- ---------- ------
   3        1  0.01396  (     zero,     0.04907)
   1     2     0.00064  (     zero,    infinity)
   1        2  0.05872  (     zero,     0.12666) **
   2     4     0.01447  (     zero,     0.04458) **
   2     5     0.00100  (     zero,    infinity)
   3     3     0.10801  (  0.01151,     0.21877) **
   3     1     0.01046  (     zero,     0.04404)

     *  = significantly positive, P < 0.05
     ** = significantly positive, P < 0.01

> plot(mltree,no.margin=TRUE,type="unrooted",edge.width=2)

OK. That's all for now. We're still working on documentation & examples, but hopefully we will have a version on CRAN before too long.

1 comment: