Ns126:Calendar/NOTES/2017-1-19
Jump to navigation
Jump to search
- R script for hapinfo2LDR2 and Dprimer and bsp-plot
library("genetics") hapinfo<-c(rep("CCCC",8),rep("TTTT",2)) # average linkage among all CpG locis in haplotype hapinfo2LDR2(hapinfo) # circle plot base on hapinfo bspplot(hapinfo2matrix(hapinfo)) # linkage of loci 1 and loci 2 diseq(genotype(metype2Geno(hapinfo)1)) metype2Geno<-function(hapinfo){ rlt<-list() for(i in 1:(nchar(hapinfo)[1]-1)){ Metype<-unlist(lapply(hapinfo,function(x) substr(x,i,i+1))) Geno<-lapply(Metype,meth2Geno) rlti<-as.vector(unlist(Geno)) } return(rlt) } hapinfo2LDR2<-function(hapinfo){ # R2 # Dpprime to indicate strength of CC and TT, not same as genetic define (CT) rlt<-list() Mhap<-lapply(metype2Geno(hapinfo),function(x){diseq(genotype(as.vector(x)))}) R2<-mean(unlist(lapply(Mhap,function(x) x$R2.overall)),na.rm=T) Dpprime<-mean(unlist(lapply(Mhap,function(x) x$Dprime.overall)),na.rm=T) rlt$R2<-R2 rlt$Dpprime<--Dpprime return(rlt) } meth2Geno<-function(x){ as.list(paste(as.character(unlist(strsplit(x,split=""))),collapse="/")) } hapinfo2matrix<-function(x){ matrix<-do.call("rbind",as.list(hapinfo)) matrix<-data.matrix(do.call("rbind",sapply(matrix,function(x) strsplit(x,"")))) matrix[matrix=="C"]=1 matrix[matrix=="T"]=0 rownames(matrix)<-make.names(rownames(matrix),unique = TRUE) matrix<-matrix(as.numeric(matrix),nrow=nrow(matrix),byrow=F) matrix } bspplot<-function(Matrix){ par(mar=c(3,3,1,5)) col=colorRampPalette(c("white", "red"))(20) circle=c(1,19) plot(x=nrow(Matrix),y=ncol(Matrix),type="n",xlab="",ylab="",xlim=c(0,ncol(Matrix)+1),ylim=c(0,nrow(Matrix)+1)) for(i in 1:ncol(Matrix)){ for(j in 1:nrow(Matrix)){ points(i,j,col=1,pch=circle[Matrix[j,i]+1],cex=1) } } }