您当前的位置:首页 > r

R语言 画曼哈顿图 修正版1.1

国庆前写了个为GWAS海量数据画曼哈顿图的函数,在后面的测试和使用中有出现了一点点问题,经过国庆的努力后稍作修改,修正了一些bug,主要是在处理数据时会时不时遇到为NA的数据,所以对此进行了一些预处理,保证程序安全和正常运行。
o22linux系统宝典
函数代码下载
o22linux系统宝典
o22linux系统宝典
代码如下:
o22linux系统宝典
o22linux系统宝典
o22linux系统宝典
o22linux系统宝典
mhtObject<-setClass("mhtObject",representation(chr="character",bp="numeric",pvalue="numeric"),prototype=prototype(chr=c(rep("chr1",100),rep("chr7",200),rep("chrX",150)),bp=c(120:219,32:231,15:164),pvalue=abs(runif(450,0,1))))
o22linux系统宝典
x<-new("mhtObject")
o22linux系统宝典
o22linux系统宝典
mhtplot<-function(data,lab.col="white",axis.col="black",bg.col="brown",col=c("white","green","blue"),tlog=FALSE,...)
{
    
    if(class(data) !="mhtObject")
      stop("It is not a S4 class /n")
     
    if(tlog == TRUE)
    {
      tmp<-log(10,data@pvalue)
      data@pvalue<--tmp
    }
    
    
   f<-levels(factor(data@chr))
   colv<-col
   col<-""
    for ( i in1:ceiling(length(f)/2))
    {
      col[2*i-1]<-rainbow(length(f),start=0.9,end=0.5)[i]
      col[2*i]<-rainbow(length(f),start=0.5,end=0.1)[ceiling(length(f)/2)-i+1]
    }
   #col<-rainbow(length(f),start=0.7,end=0.1)
    col<-rep(colv,ceiling(length(f)/length(colv)))
    maxt<-length(f)*10
    
    
    
#     plot.new()
    par(bg=bg.col,col.lab=lab.col,col.axis=axis.col,xaxs="r")
    par(lab=c(5,5,7))
    plot(rep(maxt,length(data@pvalue)),data@pvalue,type="n",xlim=c(0,maxt/10),xlab="BasePosition",ylab="P-value",axes=FALSE)
    axis(1,seq(.5,length(f)-0.5,1),labels=f,tick=FALSE)
    axis(1,seq(0,length(f),1),labels=FALSE)
    
o22linux系统宝典
#axis(2,seq(min(data@pvalue,na.rm=TRUE),max(data@pvalue,na.rm=TRUE)+(max(data@pvalue,na.rm=TRUE)-min(data@pvalue,na.rm=TRUE)*0.1),(max(data@pvalue,na.rm=TRUE)-min(x@pvalue,na.rm=TRUE))),lwd=0.5)
o22linux系统宝典
if(tlog==TRUE)
{
   axis(2,at=c(0,ceiling(max(data@pvalue,na.rm=TRUE))),lwd=0.5)
   abline(h=7,lty=3,col="white")
}
else
{
   abline(h=max(data@pvalue,na.rm=TRUE)*0.95,lty=3,col="white")
   axis(2,at=c(0,1),lwd=0.5)
o22linux系统宝典
o22linux系统宝典
}
o22linux系统宝典
o22linux系统宝典
    
o22linux系统宝典
    for ( i in 1:length(f))
    {
      k<-data@bp[data@chr==f[i]]
      n=1
      xt<-0
      for( j in 1:length(k) )
      {
        xt[n]<-(k[j]-min(k))*.8/(max(k)-min(k))+(i-1)+0.1
         n<-n+1
      }
      points(xt,data@pvalue[data@chr==f[i]],col=col[i],...)
    }
}
jpeg(width=1000,height=618)
mhtplot(x,pch=20,cex=0.8,tlog=TRUE,bg.col="light green")
o22linux系统宝典


沪ICP备10206494号-4