R语言 画曼哈顿图

A Manhattan plot is a type of scatter plot, usually used to display data with a large number of data-points – many of non-zero amplitude, and with a distribution of higher-magnitude values, for instance in genome-wide association studies(GWAS).^[1] In GWAS Manhattan plots, genomic coordinates are displayed along the X-axis, with the negative logarithm of the association P-value for each single nucleotide polymorphism displayed on the Y-axis. Because the strongest associations have the smallest P-values (e.g., 10^-15), their negative logarithms will be the greatest (e.g., 15).

It gains its name from the similarity of such a plot to the Manhattan skyline: a profile of skyscrapers towering above the lower level “buildings” which vary around a lower height

mhtplot()

使用方法：创建mhtObject S4类，该类包括三个slot，第一个为chr，为字符值，第二个为bp为碱基位置，第三个为pvalue，第二第三个均为数值型

创建命令为x<-mhtObject(chr=c(“chr1″,”chr2″,”chr3”),bp=c(1,2,3),pvalue=c(1,2,3))

mhtObject包含有一个默认的样本数据，可以通过创建空类获得： x<-mhtObject()

mhtplot()的可用参数还有：

bg.col

   背景颜色，默认为”brown”

lab.col

   坐标轴标签颜色，默认为”white”

axis.col

   坐标轴颜色，默认为”black”

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(rchisq(450,0.05,0.05))))

mhtplot<-function(data,lab.col=”white”,axis.col=”black”,bg.col=”brown”,col=c(“white”,”green”,”blue”),tlog=FALSE,…)

{

    

    if(class(data) != “mhtObject”)

       stop(“不是一个合法的mhtObject S4类”)

     

    if(tlog == TRUE)

    {

       tmp<-log(10,data@pvalue)

       data@pvalue<-tmp

    }

    

    

    f<-levels(factor(data@chr))

    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(data@pvalue)),data@pvalue,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(data@pvalue),max(data@pvalue)+(max(data@pvalue)-min(data@pvalue)*0.1),(max(data@pvalue)-min(x@pvalue))/10),lwd=0.5)

     

    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)

          n<-n+1

       }

       points(xt,data@pvalue[data@chr==f[i]],col=col[i],…)

    }

    abline(h=max(data@pvalue)*0.95,lty=3,col=”dark blue”)

}

x<-mhtObject(chr=c(rep(“chr1”,100),rep(“chr2”,200),rep(“chrX”,150)),bp=c(120:219,32:231,15:164),pvalue=abs(rchisq(450,1)))

k<-mhtplot(x,pch=20,cex=0.4)