The Manhattan map (manhattan Plot) in R language Painting Genome-wide Association analysis

1. Install R in Linux

2, ready to draw the Manhattan map of the R script is MANHATTAN.R,MANHATTAN.R content as Follows:

#!/usr/bin/rscript#example:rscript PLOT_MANHATOM.R xxx.assoc xxx.pdfargv <-commandargs () #define The function to the PLO T 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) W.H.O.Le_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)} p ng (save_file,height=500,width=600) par (mar=c (5,6,4,2)) plot (xlabel,flogp,col=whole_col,axes=f,xlab= "", ylab= "", yl Im=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") p Oints (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 AX                is character #las =0 all parallel to the axis #las =2 all horizontal for (i-1:22) {        Mtext (label[i],side=1,at=lab[i],las=0,font=1,cex=0.8) #cex Magnified relative 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 PLO T# 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&LT;-ARGV[6];p rint (site_file) save_file<-argv[7];p rint (save_file) plot_manhatton (site_file,save_file) #round (x,digits=3 ) Keep the length of the digit

3, ready to Plink run out Of. assoc.linear files, such as Mydata.assoc.linear

4. Enter one of the following commands in Linux:

Among them, Mydata.png is the Manhattan map we want (manhattan Plot)

Rscript MANHATTAN.R mydata.assoc.linear Mydata.png

The Manhattan map (manhattan Plot) in R language Painting Genome-wide Association analysis