R语言画全基因组关联分析中的曼哈顿图（manhattan plot）

#!/usr/bin/Rscript
#example : Rscript plot_manhatom.r XXX.assoc XXX.pdf
argv <- commandArgs()
#define the function to plot the manhatton and quantitle-quantitle plot 
plot_manhatton<-function(site_file,save_file){
        read.table(site_file,header=T)->data
    data<-data[which(data[,5]=="ADD"),]
    logp=-log10(data[,9]) #value
        chr=data[,1] #chr
        position<-data[,3] #position
#       xaxis=0
        sig=numeric()
        lab=numeric()
    flogp=numeric();
        the_col=c("darkblue","gray"); #chr class
        whole_col=c()
    xlabel=numeric();
        length_add=0
        label=numeric();
        for(i in 1:22){
                position[chr==i]->chr_pos
        chr_pos<-chr_pos-chr_pos[1]+17000000
        chr_num=length(chr_pos)
                cat("For chrosome",i,",Num of SNPs:",chr_num,"\n")
                if(length(chr_pos)==0){
                        next
                }
                flogp=c(flogp,logp[chr==i])
                label=c(label,i)
                whole_col=c(whole_col,rep(the_col[i%%2+1],chr_num))
                chr_pos=length_add+chr_pos
                xlabel=c(xlabel,chr_pos)
                length_add=sort(chr_pos,decreasing=T)[1]
        lab=c(lab,(chr_pos[1]+length_add)/2)
        }
    png(save_file,height=500,width=600)
    par(mar=c(5,6,4,2))
        plot(xlabel,flogp,col=whole_col,axes=F,xlab="",ylab="",ylim=c(0,12),pch=20,cex=0.5,cex.lab=1.2,cex.axis=1.4)#ylim=c(0,6)
    significance=-log10(0.05/length(data[,1]))
    sig_pos<-xlabel[which(flogp>significance)]
    for(i in 1:length(sig_pos)){
        sig<-c(sig,which(xlabel>sig_pos[i]-60000 & xlabel<sig_pos[i]+60000))
    }
    sig<-unique(sig)
    cat("significant signal:",sig[10:20],"\n")
    points(xlabel[sig],flogp[sig],col="red",pch=20,cex=0.5)
    title(xlab="Chromosome",ylab=expression(-log[10]*(p)),cex.lab=1.4)
#   title(xlab="Chromosome",ylab="Score",cex.lab=1.8)
#       title("GWAS scan for Hair curl", cex.main=2.5)
#       yaxis=seq(0,1,by=0.1)
#       axis(2,yaxis,yaxis,las=1,cex.axis=1.5,line=-2)
        axis(2,las=2,cex.axis=1.4)
        #las can change the directory of the axis character
        #las=0 always parallel to the axis
        #las=2 always horizontal


        for(i in 1:22){
                mtext(label[i],side=1,at=lab[i],las=0,font=1,cex=0.8)
                #cex magnified relative to the to the default
        }
#    mtext("X",side=1,at=x[23],las=0,font=1,cex=1.4)
#    mtext("Y",side=1,at=x[24],las=0,font=1,cex=1.4)

#       axis(1,x,as.character(label),tick=TRUE,font=1)
        par(lty=2)
        abline(h=significance,cex=1.5,col="red")
#plot the QQ plot
#    par(fig=c(0.58,0.92,0.43,0.95),new=TRUE)
#    observed <- sort(data[,12])
#    logp=-log10(observed) 
#    expected <- c(1:length(observed))
#    lexp <- -(log10(expected / (length(expected)+1)))
#    plot(x=lexp,y=logp,pch=19,cex=0.6, xlab="Expected (-logP)", ylab="Observed (-logP)",cex.lab=1.2,cex.axis=1.2)
#    abline(0,1,col="red",lwd=2,lty=1)

    dev.off()
}


site_file<-argv[6];print(site_file)
save_file<-argv[7];print(save_file)
plot_manhatton(site_file,save_file)


#round(x,digits=3) keep the length of the digit

Rscript manhattan.r mydata.assoc.linear mydata.png