RG#27: Manhattan plot

# data 

set.seed(123)
data22 <- data.frame (snp = 1: 20000*5 , chr = c(rep(1:5, each = 20000)), 
pos= rep(1:20000, 5), pval1= rnorm(20000*5, 0.2, 0.3), 
pval2 = rnorm(20000*5, 0.2, 0.3) )
# the above simulation produce negative values so the following will replace 
# negative values with NA
data22$pval1[data22$pval1 < 0] <- NA
# removal of negative values 
dat2 <- data22[!is.na(data22$pval1),]

# install plantbreeding program from Rforge

require(plantbreeding)   
op <- par(mfrow=c(2,1), cex.axis = 0.75, font = 1, family = "serif")
par(op)

manhatton.plot(dataframe = dat2, SNPname = "snp", chromosome = "chr", position = "pos", pvcol = "pval1", line1 = 4, line2 = 8, ymax = "maximum", ymin = 0, gapbp = 2000)

title(main = "Mahattan plot of results for trait1", sub = "Method: Linear mixed model")

manhatton.plot(dataframe = dat2, SNPname = "snp", chromosome = "chr", position = "pos", pvcol = "pval2", line1 = 4, line2 = 8, ymax = "maximum", ymin = 0, color = c("gray50", "gray20"), pch = c(19, 1))

title(main = "Mahattan plot of results for trait2", sub = "Method: Linear mixed model")

data12 <- data.frame (snp = 1: 2000*20 , chr = c(rep(1:20, each = 2000)), 
pos= rep(1:2000, 20), pval= rnorm(2000*20, 0.001, 0.005))

# adding stripes (color filled rectangles)
manhatton.plot(dataframe = data12, SNPname = "snp", chromosome = "chr", 
position = "pos", pvcol = "pval",ymax = "maximum", ymin = 0,  gapbp = 500, 
color=c("hotpink3","dodgerblue4"), line1 = 3, line2 = 5, pch = c(1,20),
stripe= TRUE, stripe.color = c("aliceblue", "cornsilk", "lightgoldenrod") )

# with annotation 
manhatton.plot(dataframe = data12, SNPname = "snp", chromosome = "chr", 
position = "pos", pvcol = "pval",ymax = "maximum", ymin = 0,  gapbp = 500, 
color=c("hotpink3","dodgerblue4"), line1 = 3, line2 = 5, pch = c(1,20),
annotate= TRUE, threshold= 6, thresholddir = "upper" )

# plot R2 or p-value 

# Example 3, plotting simple p-value or R-square 
 data13 <- data.frame (snp = 1: 2000*20 , chr = c(rep(1:20, each = 2000)), 
pos= rep(1:2000, 20), pval= rnorm(2000*20, 0.15, 0.05))
manhatton.plot(dataframe = data13, pconv= "NULL", ylabel = "R-square", 
SNPname = "snp", chromosome = "chr", 
position = "pos", pvcol = "pval",ymax = "maximum", ymin = 0,  gapbp = 500, 
color=c("hotpink3","dodgerblue4"), line1 = 0.6, line2 = 0.2, pch = c(1,20) )