Simple DNA sequence simulator using sim.history internally

I was recently playing with simSeq in the phangorn package - but I couldn't get it to do exactly what I wanted (probably for lack of sufficient patience). Then I realized that I could (nearly) just as easily simulate DNA sequences using phytools with the function sim.history. Here's a quick & incredibly simple function that I wrote to do this that wraps around sim.history:

genSeq<-function(tree,l=1000,Q=NULL,rate=1, format="DNAbin",...){
  if(is.null(Q)){
    Q<-matrix(1,4,4)
    rownames(Q)<-colnames(Q)<-c("a","c","g","t")
    diag(Q)<-0
    diag(Q)<--colSums(Q)
  }
  X<-replicate(l,sim.history(tree,rate*Q)$states)
  if(format=="DNAbin") return(as.DNAbin(X))
  else if(format=="phyDat") return(as.phyDat(X))
  else if(format=="matrix") return(X)
}

Yes, it's really that simple. Now, admittedly this function cannot simulate rate heterogeneity among sites, unequal base frequencies, or invariant sites, but these would be relatively straightforward to add.

Here's, let's try it out:

> ## first let's generate a random tree
> ## with a basal root taxon "Z"
> tree<-pbtree(n=25,scale=0.9)
> tree$tip.label<-LETTERS[25:1]
> tree$root.edge<-0
> root<-list(edge=matrix(c(3,1,3,2),2,2,byrow=TRUE),
edge.length=c(1,0.1),
tip.label=c("Z","NA"),
Nnode=1)
> class(root)<-"phylo"
> tree<-paste.tree(root,tree)
>
> ## now let's simulate under Juke-Cantor
> ## (the default)
> X<-genSeq(tree,l=2000,rate=0.1,format="phyDat")
> X
26 sequences with 2000 character and 1711 different site patterns.
The states are a c g t
>
> ## now let's use our data for inference
> library(phangorn)
> obj<-pml(rtree(n=26,tip.label=LETTERS),X)
> fit<-optim.pml(obj,optNni=TRUE)
optimize edge weights: -56156.91 --> -40367.43
optimize edge weights: -40367.43 --> -40367.43
optimize topology: -40367.43 --> -38967.95
...
optimize topology: -25858.87 --> -25858.87
0
Warning message:
I unrooted the tree (rooted trees are not yet supported)
>
> ## measure RF distance to original
> RF.dist(unroot(tree),unroot(fit$tree))
[1] 0
> ## plot original and estimated tree
> par(mfcol=c(2,1))
> plotTree(tree,mar=c(0.1,1.1,3.1,0.1))
> title("a) Generating tree",adj=0,cex.main=1.2)
> plotTree(midpoint(fit$tree),mar=c(0.1,1.1,3.1,0.1))
> title("b) Estimated tree using ML",adj=0,cex.main=1.2)

That works OK. Adding other attributes typical of molecular sequence models is really easy, so I'll probably do that. Pretty cool.