nievergeltlab
3/13/2017 - 4:12 PM
Manhattan Plots V2 of the plotter will color SNPs within a set kb window of any hits
Manhattan Plots
V2 of the plotter will color SNPs within a set kb window of any hits
#############################Begin Function#############################
#First we define the function with required parameters and optional parameters with their defaults
ManhattanPlot_AJS_cut <- function(chr, pos, pvalue, snp, genomesig = 0, genomesug = 0,photoname = "Manhattan Plot", outname = "ManhattanPlot_AJS", colors = c("navyblue","gray69"), sigcolor = "darkred", sugcolor = "indianred", ncex = .5, plotsymbol = 16, sugcex = 1, sigcex = 1.5, sigsnpcolor = "red", nonautosome = c(23,24,25,26),xlabel = "Chromosomal Positions",ylabel = "-log10(p-value)", pdf = "TRUE",topsnplabels = "FALSE", pvalue_miss = "NA") {
#bind the input data into one table and prune out missing data based on the pvalue_miss option and then
#redefines the input variables chr, pos, pvalue, snp to reflect the cleaned data
data_cl <- data.frame(snp,chr,pos,pvalue)
data_cl <- subset(data_cl,data_cl$pvalue != pvalue_miss)
chr <- data_cl$chr
pos <- data_cl$pos
pvalue <- data_cl$pvalue
snp <- data_cl$snp
#convert genomesig and genomsug pvalue options into the –log10 form for use on the plot
if (genomesug != 0) {
genomesug <- -log10(genomesug)
}
if (genomesig != 0) {
genomesig <- -log10(genomesig)
}
#create a list of unique chromosome codes
chroms <- as.numeric(levels(as.factor(chr)))
#create a table with important chromosome data for each unique chromosome code.
chrominfo <- data.frame()
colflag <- 0
#loop through each unique chromosome code and store in a table the code [1], the max bp position [2],
#the cumulative bp position to mark the beginning of the chromosome on the plot [3], the cumulative
#midpoint of each chromosome to draw the label at [4], the label to be used on the plot for each
#chromosome [5], the color code for the chromosome [6]
for (i in 1:length(chroms)) {
chrominfo[i,1] <- chroms[i]
chrominfo[i,2] <- max(pos[chr == chroms[i]], na.rm = T)
ifelse(i == 1, chrominfo[i,3] <- 0, chrominfo[i,3] <- chrominfo[i-1,2] + chrominfo[i-1,3])
chrominfo[i,4] <- chrominfo[i,3] + (.5*chrominfo[i,2])
if (chroms[i] <= 22) {
chrominfo[i,5] <- chroms[i]
} else if (chroms[i] == nonautosome[1]) {
chrominfo[i,5] <- "X"
} else if (chroms[i] == nonautosome[2]) {
chrominfo[i,5] <- "Y"
} else if (chroms[i] == nonautosome[3]) {
chrominfo[i,5] <- "XY"
} else if (chroms[i] == nonautosome[4]) {
chrominfo[i,5] <- "MT"
}
if (colflag == 0) {
chrominfo[i,6] <- colors[1]
colflag <- 1
} else {
chrominfo[i,6] <- colors[2]
colflag <- 0
}
#5/2/13 AM edit, this makes it so that the unplaced stuff is ALWAYS grey, so the color choice is uniform (meta analysis files do not always have unplaced SNPs)
if (chroms[i] == 0)
{
chrominfo[i,6] <- colors[2]
colflag <- 0
}
}
#define a list of the sizes of each point on the graph based on its significance level and the
#defined cex for ncex, sugcex, and sigcex and the defined genomesig and genomesug options
ptcex <- pvalue
for (i in 1:length(pvalue)){
if (genomesig == 0) {
if (genomesug == 0) {
ptcex <- ncex
}else {
if (-log10(pvalue[i]) > genomesug) {
ptcex[i] <- sugcex
}else (ptcex[i] <- ncex)
}
}else {
if (-log10(pvalue[i]) > genomesig) {
ptcex[i] <- sigcex
}else {
if (genomesug == 0) {
ptcex[i] <- ncex
}else {
if (-log10(pvalue[i]) > genomesug) {
ptcex[i] <- sugcex
} else (ptcex[i] <- ncex)
}
}
}
}
#create an output file, either pdf or jpeg based on specified pdf option
if (pdf == "TRUE") {
pdf(file = paste(outname, ".pdf", sep = ""),11,8.5)
} else {
jpeg(filename = paste(outname, ".jpeg", sep = ""),width = 11, height = 8.5, units = "in", res = 500,quality = 100)
}
#define the maximum value along the x chromosome to be used for defining plot and axis sizes
xmax <- chrominfo[length(chrominfo[,2]),2] + chrominfo[length(chrominfo[,3]),3]
#defining the plot window size within the output file by specifying margin sizes and mute automatic
#labels to allow custom labels
#bty set to l to remove boundaries on upper right
par(mar=c(5,5,5,5),xaxs = "i", yaxs = "i", bty="l", las=1)
#Adam: Use the default PCH for all markers execept the significant ones. For these, use a bordered circle (color set by sigsnpcolor in parameters). the bg= part of the plot command will take care of the color param
#identify proximal markers, color those too
sigmarkers <- which(pvalue <= 5e-8) #Adam
posfit <- rep(FALSE,length(pos)) #Turn to true if its NEAR a top hit
boundary_bp=50000
for (i in c(1:length(pos)))
{
c1 <- pos[i] >= pos[sigmarkers] - boundary_bp
c2 <- pos[i] <= pos[sigmarkers] + boundary_bp
c3 <- chr[i] == chr[sigmarkers]
summed <- c1 + c2 + c3
if(max(summed) == 3)
{
posfit[i] <- TRUE
}
}
sigmarkers <- which(posfit == TRUE) #can just use the same variable, since we test each marker, they will be included in this DF.
marker_colors <- chrominfo[,6][match(chr,chrominfo[,1])]
marker_colors[sigmarkers] <- "black"
marker_pch <- rep(plotsymbol,length(pos))
marker_pch[sigmarkers] <- 21
#plot the points and define size of axes
plot(pos + chrominfo[,3][match(chr,chrominfo[,1])],-log10(pvalue), bg=sigsnpcolor, col = marker_colors,cex = ptcex, pch = marker_pch, ann = F, xaxt = "n", xlim = c(0,1.015*xmax), ylim = c(2, max(10)),cex.axis=1.45, cex.lab=1.5)
#Replot, where the non significant markers are NOT plotted
marker_colors[which(posfit != TRUE)] <- NA
points(pos + chrominfo[,3][match(chr,chrominfo[,1])],-log10(pvalue), bg=sigsnpcolor, col = marker_colors,cex = ptcex, pch = marker_pch)
#label the plot and axes
title(main = photoname, col.main = colors[1], cex.main = 2, font.main = 2, xlab = xlabel,ylab = ylabel,cex.lab=1.5, cex.axis=1.45 )
#draw and label axis tick marks on x axis
axis(1, at = chrominfo[,4],labels = chrominfo[,5],cex.axis=1.3, cex.lab=1.5)
#draw significant and suggestive lines according to genomesig and genomesug
genomesigcheck <- 1
genomesugcheck <- 1
if (genomesig != 0) {
if (TRUE | (max(-log10(pvalue)) > genomesig)) {
abline(h=genomesig,lty = 1, col = sigcolor)
text(1.0075*xmax,genomesig,label = "",col = sigcolor, xpd = T, pos = 4)
genomesigcheck <- 1
}
}
if (genomesug != 0) {
if (TRUE | (max(-log10(pvalue)) > genomesug)) {
abline(h=genomesug,lty = 2, col = sugcolor) #Adam edit: no suggestive line
text(1.0075*xmax,genomesug,label = "",col = sugcolor, xpd = T, pos = 4)
genomesugcheck <- 1
}
}
#label snps according to topsnplabels
if (topsnplabels == "ALL") {
for (i in chrominfo[,1]) {
textxcord <- pos[pvalue == min(pvalue[chr == i])] + chrominfo[,3][chrominfo[,1] == i]
textycord <- -log10(min(pvalue[chr == i]))
textlabel <- snp[pvalue == min(pvalue[chr == i])]
textpvalue <- min(pvalue[chr == i])
if (genomesigcheck == 1) {
if (textycord > genomesig) {
text(textxcord,textycord,textlabel,pos = 3,cex = sigcex*(2/3), col = sigcolor, offset = sigcex, xpd = T)
text(textxcord,textycord,paste("p=",textpvalue),pos = 3,cex = sigcex*(2/3), col = sigcolor, xpd = T)
}
}
if (genomesugcheck == 1) {
if (genomesig == 0) {
if (textycord > genomesug) {
text(textxcord,textycord,textlabel,pos = 3,cex = sugcex*(2/3), col = sugcolor, offset = sugcex, xpd = T)
text(textxcord,textycord,paste("p=",textpvalue),pos = 3,cex = sugcex*(2/3), col = sugcolor, xpd = T)
}
}
else {
if (textycord > genomesug) {
if (textycord < genomesig) {
text(textxcord,textycord,textlabel,pos = 3,cex = sugcex*(2/3), col = sugcolor, offset = sugcex, xpd = T)
text(textxcord,textycord,paste("p=",textpvalue),pos = 3,cex = sugcex*(2/3), col = sugcolor, xpd = T)
}
}
}
}
}
}
else {
if (topsnplabels == "TOP") {
textxcord <- pos[pvalue == min(pvalue)] + chrominfo[,3][chrominfo[,1] == chr[pvalue == min(pvalue)]]
textycord <- -log10(min(pvalue))
textlabel <- snp[pvalue == min(pvalue)]
textpvalue <- min(pvalue)
if (genomesigcheck == 1) {
if (textycord > genomesig) {
text(textxcord,textycord,textlabel,pos = 3,cex = sigcex*(2/3), col = sigcolor, offset = sigcex, xpd = T)
text(textxcord,textycord,paste("p=",textpvalue),pos = 3,cex = sigcex*(2/3), col = sigcolor, xpd = T)
}
}
if (genomesugcheck == 1) {
if (genomesig == 0) {
if (textycord > genomesug) {
text(textxcord,textycord,textlabel,pos = 3,cex = sugcex*(2/3), col = sugcolor, offset = sugcex, xpd = T)
text(textxcord,textycord,paste("p=",textpvalue),pos = 3,cex = sugcex*(2/3), col = sugcolor, xpd = T)
}
}
else {
if (textycord > genomesug) {
if (textycord < genomesig) {
text(textxcord,textycord,textlabel,pos = 3,cex = sugcex*(2/3), col = sugcolor, offset = sugcex, xpd = T)
text(textxcord,textycord,paste("p=",textpvalue),pos = 3,cex = sugcex*(2/3), col = sugcolor, xpd = T)
}
}
}
}
}
else {
if (topsnplabels == "SIG") {
for (i in chrominfo[,1]) {
textxcord <- pos[pvalue == min(pvalue[chr == i])] + chrominfo[,3][chrominfo[,1] == i]
textycord <- -log10(min(pvalue[chr == i]))
textlabel <- snp[pvalue == min(pvalue[chr == i])]
textpvalue <- min(pvalue[chr == i])
if (genomesigcheck == 1) {
if (textycord > genomesig) {
text(textxcord,textycord,textlabel,pos = 3,cex = sigcex*(2/3), col = sigcolor, offset = sigcex, xpd = T)
text(textxcord,textycord,paste("p=",textpvalue),pos = 3,cex = sigcex*(2/3), col = sigcolor, xpd = T)
}
}
}
}
}
}
#close output file
dev.off()
}
#############################End Function#############################
#############################Begin Function#############################
#First we define the function with required parameters and optional parameters with their defaults
ManhattanPlot_AJS_cut <- function(chr, pos, pvalue, snp, genomesig = 0, genomesug = 0,photoname = "Manhattan Plot", outname = "ManhattanPlot_AJS", colors = c("navyblue","gray69"), sigcolor = "darkred", sugcolor = "indianred", ncex = .5, plotsymbol = 18, sugcex = 1, sigcex = 1.5, nonautosome = c(23,24,25,26),xlabel = "Chromosomal Positions",ylabel = "-log10(p-value)", pdf = "TRUE",topsnplabels = "FALSE", pvalue_miss = "NA") {
#bind the input data into one table and prune out missing data based on the pvalue_miss option and then
#redefines the input variables chr, pos, pvalue, snp to reflect the cleaned data
data_cl <- data.frame(snp,chr,pos,pvalue)
data_cl <- subset(data_cl,data_cl$pvalue != pvalue_miss)
chr <- data_cl$chr
pos <- data_cl$pos
pvalue <- data_cl$pvalue
snp <- data_cl$snp
#convert genomesig and genomsug pvalue options into the –log10 form for use on the plot
if (genomesug != 0) {
genomesug <- -log10(genomesug)
}
if (genomesig != 0) {
genomesig <- -log10(genomesig)
}
#create a list of unique chromosome codes
chroms <- as.numeric(levels(as.factor(chr)))
#create a table with important chromosome data for each unique chromosome code.
chrominfo <- data.frame()
colflag <- 0
#loop through each unique chromosome code and store in a table the code [1], the max bp position [2],
#the cumulative bp position to mark the beginning of the chromosome on the plot [3], the cumulative
#midpoint of each chromosome to draw the label at [4], the label to be used on the plot for each
#chromosome [5], the color code for the chromosome [6]
for (i in 1:length(chroms)) {
chrominfo[i,1] <- chroms[i]
chrominfo[i,2] <- max(pos[chr == chroms[i]], na.rm = T)
ifelse(i == 1, chrominfo[i,3] <- 0, chrominfo[i,3] <- chrominfo[i-1,2] + chrominfo[i-1,3])
chrominfo[i,4] <- chrominfo[i,3] + (.5*chrominfo[i,2])
if (chroms[i] <= 22) {
chrominfo[i,5] <- chroms[i]
} else if (chroms[i] == nonautosome[1]) {
chrominfo[i,5] <- "X"
} else if (chroms[i] == nonautosome[2]) {
chrominfo[i,5] <- "Y"
} else if (chroms[i] == nonautosome[3]) {
chrominfo[i,5] <- "XY"
} else if (chroms[i] == nonautosome[4]) {
chrominfo[i,5] <- "MT"
}
if (colflag == 0) {
chrominfo[i,6] <- colors[1]
colflag <- 1
} else {
chrominfo[i,6] <- colors[2]
colflag <- 0
}
#5/2/13 AM edit, this makes it so that the unplaced stuff is ALWAYS grey, so the color choice is uniform (meta analysis files do not always have unplaced SNPs)
if (chroms[i] == 0)
{
chrominfo[i,6] <- colors[2]
colflag <- 0
}
}
#define a list of the sizes of each point on the graph based on its significance level and the
#defined cex for ncex, sugcex, and sigcex and the defined genomesig and genomesug options
ptcex <- pvalue
for (i in 1:length(pvalue)){
if (genomesig == 0) {
if (genomesug == 0) {
ptcex <- ncex
}else {
if (-log10(pvalue[i]) > genomesug) {
ptcex[i] <- sugcex
}else (ptcex[i] <- ncex)
}
}else {
if (-log10(pvalue[i]) > genomesig) {
ptcex[i] <- sigcex
}else {
if (genomesug == 0) {
ptcex[i] <- ncex
}else {
if (-log10(pvalue[i]) > genomesug) {
ptcex[i] <- sugcex
} else (ptcex[i] <- ncex)
}
}
}
}
#create an output file, either pdf or jpeg based on specified pdf option
if (pdf == "TRUE") {
pdf(file = paste(outname, ".pdf", sep = ""),11,8.5)
} else {
jpeg(filename = paste(outname, ".jpeg", sep = ""),width = 11, height = 8.5, units = "in", res = 500,quality = 100)
}
#define the maximum value along the x chromosome to be used for defining plot and axis sizes
xmax <- chrominfo[length(chrominfo[,2]),2] + chrominfo[length(chrominfo[,3]),3]
#defining the plot window size within the output file by specifying margin sizes and mute automatic
#labels to allow custom labels
#bty set to l to remove boundaries on upper right
par(mar=c(5,5,5,5),xaxs = "i", yaxs = "i", bty="l", las=1)
#plot the points and define size of axes
plot(pos + chrominfo[,3][match(chr,chrominfo[,1])],-log10(pvalue), col = chrominfo[,6][match(chr,chrominfo[,1])],cex = ptcex, pch = plotsymbol, ann = F, xaxt = "n", xlim = c(0,1.015*xmax), ylim = c(0, max(10.5)),cex.axis=1.45, cex.lab=1.5)
#label the plot and axes
title(main = photoname, col.main = colors[1], cex.main = 2, font.main = 2, xlab = xlabel,ylab = ylabel,cex.lab=1.5, cex.axis=1.45 )
#draw and label axis tick marks on x axis
axis(1, at = chrominfo[,4],labels = chrominfo[,5],cex.axis=1.3, cex.lab=1.5)
#draw significant and suggestive lines according to genomesig and genomesug
genomesigcheck <- 1
genomesugcheck <- 1
if (genomesig != 0) {
if (TRUE | (max(-log10(pvalue)) > genomesig)) {
abline(h=genomesig,lty = 1, col = sigcolor)
text(1.0075*xmax,genomesig,label = "",col = sigcolor, xpd = T, pos = 4)
genomesigcheck <- 1
}
}
if (genomesug != 0) {
if (TRUE | (max(-log10(pvalue)) > genomesug)) {
abline(h=genomesug,lty = 2, col = sugcolor)
text(1.0075*xmax,genomesug,label = "",col = sugcolor, xpd = T, pos = 4)
genomesugcheck <- 1
}
}
#label snps according to topsnplabels
if (topsnplabels == "ALL") {
for (i in chrominfo[,1]) {
textxcord <- pos[pvalue == min(pvalue[chr == i])] + chrominfo[,3][chrominfo[,1] == i]
textycord <- -log10(min(pvalue[chr == i]))
textlabel <- snp[pvalue == min(pvalue[chr == i])]
textpvalue <- min(pvalue[chr == i])
if (genomesigcheck == 1) {
if (textycord > genomesig) {
text(textxcord,textycord,textlabel,pos = 3,cex = sigcex*(2/3), col = sigcolor, offset = sigcex, xpd = T)
text(textxcord,textycord,paste("p=",textpvalue),pos = 3,cex = sigcex*(2/3), col = sigcolor, xpd = T)
}
}
if (genomesugcheck == 1) {
if (genomesig == 0) {
if (textycord > genomesug) {
text(textxcord,textycord,textlabel,pos = 3,cex = sugcex*(2/3), col = sugcolor, offset = sugcex, xpd = T)
text(textxcord,textycord,paste("p=",textpvalue),pos = 3,cex = sugcex*(2/3), col = sugcolor, xpd = T)
}
}
else {
if (textycord > genomesug) {
if (textycord < genomesig) {
text(textxcord,textycord,textlabel,pos = 3,cex = sugcex*(2/3), col = sugcolor, offset = sugcex, xpd = T)
text(textxcord,textycord,paste("p=",textpvalue),pos = 3,cex = sugcex*(2/3), col = sugcolor, xpd = T)
}
}
}
}
}
}
else {
if (topsnplabels == "TOP") {
textxcord <- pos[pvalue == min(pvalue)] + chrominfo[,3][chrominfo[,1] == chr[pvalue == min(pvalue)]]
textycord <- -log10(min(pvalue))
textlabel <- snp[pvalue == min(pvalue)]
textpvalue <- min(pvalue)
if (genomesigcheck == 1) {
if (textycord > genomesig) {
text(textxcord,textycord,textlabel,pos = 3,cex = sigcex*(2/3), col = sigcolor, offset = sigcex, xpd = T)
text(textxcord,textycord,paste("p=",textpvalue),pos = 3,cex = sigcex*(2/3), col = sigcolor, xpd = T)
}
}
if (genomesugcheck == 1) {
if (genomesig == 0) {
if (textycord > genomesug) {
text(textxcord,textycord,textlabel,pos = 3,cex = sugcex*(2/3), col = sugcolor, offset = sugcex, xpd = T)
text(textxcord,textycord,paste("p=",textpvalue),pos = 3,cex = sugcex*(2/3), col = sugcolor, xpd = T)
}
}
else {
if (textycord > genomesug) {
if (textycord < genomesig) {
text(textxcord,textycord,textlabel,pos = 3,cex = sugcex*(2/3), col = sugcolor, offset = sugcex, xpd = T)
text(textxcord,textycord,paste("p=",textpvalue),pos = 3,cex = sugcex*(2/3), col = sugcolor, xpd = T)
}
}
}
}
}
else {
if (topsnplabels == "SIG") {
for (i in chrominfo[,1]) {
textxcord <- pos[pvalue == min(pvalue[chr == i])] + chrominfo[,3][chrominfo[,1] == i]
textycord <- -log10(min(pvalue[chr == i]))
textlabel <- snp[pvalue == min(pvalue[chr == i])]
textpvalue <- min(pvalue[chr == i])
if (genomesigcheck == 1) {
if (textycord > genomesig) {
text(textxcord,textycord,textlabel,pos = 3,cex = sigcex*(2/3), col = sigcolor, offset = sigcex, xpd = T)
text(textxcord,textycord,paste("p=",textpvalue),pos = 3,cex = sigcex*(2/3), col = sigcolor, xpd = T)
}
}
}
}
}
}
#close output file
dev.off()
}
#############################End Function#############################args <- commandArgs(trailingOnly = TRUE)
scriptloc <- args[1]
results <- args[2]
outfile <- args[3]
snpcol <- args[4]
chrcol <- args[5]
bpcol <- args[6]
pcol <- args[7]
print ("Reading data")
#Read in data
#get number of rows in file (reads file faster)
nr <- as.numeric(system(paste("wc -l", results, " | awk \'{print $1}\' "),intern=T))
#get list of numeric column classes (reads file faster)
dat_temp <- read.table(results, header=T,nrows=50,stringsAsFactors=F)
classes_to_use <- sapply(dat_temp , class)
#read file
dat <- read.table (results, header=T,stringsAsFactors=F)
print ("Loading Manhattan Script")
source(scriptloc)
SNP <- dat[,snpcol]
CHR <- dat[,chrcol]
BP <- dat[,bpcol]
P <- dat[,pcol]
print ("Plotting data")
ManhattanPlot_AJS_cut(CHR, BP, P, SNP, genomesig = 5*(10^-8), genomesug = 5*(10^-6),photoname = '',
outname = outfile, colors = c(rgb(128,161,103,max=255),'black'), sigcolor = 'darkred', sugcolor = 'indianred', ncex = .5,
plotsymbol = 18, sugcex = 1, sigcex = 1.5, nonautosome = c(23,24,25,26),xlabel = 'Chromosomal Positions',ylabel = '-log10(p-value)',
pdf = 'FALSE',topsnplabels = 'FALSE', pvalue_miss = 'NA')