In other R news, there's a cool site, R-bloggers, that is a portal to a number of other blogs that deal with R. It's great to see what other people manage to do in R and a good way to learn about its capabilities.
Happy New Year!
library(ape)
#Input: dat---an object of class 'DNAbin'
titv<-function(dat){
mat<-as.matrix(dat)
res<-matrix(NA, ncol=dim(mat)[1], nrow=dim(mat)[1], dimnames=list(x=names(dat), y=names(dat)))
for(i in 1:dim(mat)[1]){
for(j in 1:dim(mat)[1]){
vec<-as.numeric(mat[i,])+as.numeric(mat[j,])-8
res[i,j]<-length(grep("200|56",vec)) #Transitions
res[j,i]<-length(grep("152|168|88|104",vec)) #Transversions
}
}
res
}
#Example
data(woodmouse)
ti<-titv(woodmouse)
tv<-t(ti)
tv[lower.tri(tv)] #Number of transversions
ti[lower.tri(ti)] #Number of transitions
#Saturation plot
dist<-dist.dna(woodmouse)
plot(ti[lower.tri(ti)]~dist)
points(tv[lower.tri(tv)]~dist, pch=20, col="red")
2 comments:
Hi there,
Thank you for the post, I would like to suggest two things:
1) Please consider adding the link to R Bloggers:
http://www.r-bloggers.com/
2) Would you like to add your blog (the posts tagged with R) to R-bloggers?
If so, you could e-mail me about it or filling the form on:
http://www.r-bloggers.com/add-your-blog/
Best :)
Tal
This version of titv() was found to have a minor bug in it that considered the difference between an N and a T to be a transition. A faster version of the functions that doesn't give this result has been generously given to me by Klaus Schliep and can be found here.
Post a Comment