R语言 画曼哈顿图

国庆前写了个为GWAS海量数据画曼哈顿图的函数，在后头的测试和利用中有呈现了一点点问题，颠末国庆的尽力后稍作修改，批改了一些bug，主要是在处理惩罚数据时会时不时碰着为NA的数据，所以对此举办了一些预处理惩罚，担保措施安详和正常运行。

函数代码下载

代码如下：

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))))

x<-new(“mhtObject”)

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,[email protected])
[email protected]<–tmp
}


f<-levels(factor([email protected]))
colv<-col
col<-“”
for ( i in 1: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([email protected])),[email protected],type=”n”,xlim=c(0,maxt/10),xlab=”Base
Position”,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)


#axis(2,seq(min([email protected],na.rm=TRUE),max([email protected],na.rm=TRUE)+(max([email protected],na.rm=TRUE)-min([email protected],na.rm=TRUE)*0.1),(max([email protected],na.rm=TRUE)-min([email protected],na.rm=TRUE))),lwd=0.5)

if(tlog==TRUE)
{

axis(2,at=c(0,ceiling(max([email protected],na.rm=TRUE))),lwd=0.5)
abline(h=7,lty=3,col=”white”)
}
else
{

abline(h=max([email protected],na.rm=TRUE)*0.95,lty=3,col=”white”)
axis(2,at=c(0,1),lwd=0.5)

}

for ( i in 1:length(f) )
{
k<[email protected][[email protected]==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,[email protected][[email protected]==f[i]],col=col[i],…)
}
}
jpeg(width=1000,height=618)
mhtplot(x,pch=20,cex=0.8,tlog=TRUE,bg.col=”light green”)