Mflick at 14:49, 14 January 2014

2014-01-14T14:49:36Z

← Older revision		Revision as of 10:49, 14 January 2014
Line 416:		Line 416:
	should work well.		should work well.

	~~= Alternatives =~~

	Note that Cristen Willer and Serena Sanna have also made many beautiful plots for GWAS data using R base graphics. They have provided a sample script that does many of the same things as the above function, plus takes care of reading in the data, highlighting known SNPs, and making QQ plots. The sample code can be found at
	~~snowwhite:/home/cristen/scripts/makeplot_manhattan_qq_forsharing.R~~

	[[Category:Code Samples]]		[[Category:Code Samples]]

Goncalo at 19:53, 9 January 2013

2013-01-09T19:53:17Z

← Older revision		Revision as of 15:53, 9 January 2013
Line 2:		Line 2:

	[[File:Sample-manhattan-plot.png\|center]]		[[File:Sample-manhattan-plot.png\|center]]

			= Credit =

			This page is based on a tutorial originally written by [mailto:mflick@umich.edu Matthew Flickinger].

	= Manhattan Plotting Function =		= Manhattan Plotting Function =

Goncalo: Created page with 'A useful way to summarize genome-wide association data is with a Manhattan plot. This type of plot has a point for every SNP or location tested with the position in the genome al…'

2013-01-09T19:17:55Z

Created page with 'A useful way to summarize genome-wide association data is with a Manhattan plot. This type of plot has a point for every SNP or location tested with the position in the genome al…'

New page

A useful way to summarize genome-wide association data is with a Manhattan plot. This type of plot has a point for every SNP or location tested with the position in the genome along the x-axis and the -log10 p-value on the y-axis.

[[File:Sample-manhattan-plot.png|center]]

= Manhattan Plotting Function =

Here is a function which can make a Manhattan plot using lattice graphics. While the function itself is quite long, you don't have to worry about most of it. You really only need to pay attention to the parameters that you pass to the funciton. An example of its use is given below.
<syntaxhighlight lang="rsplus">
library(lattice)
manhattan.plot<-function(chr, pos, pvalue,
sig.level=NA, annotate=NULL, ann.default=list(),
should.thin=T, thin.pos.places=2, thin.logp.places=2,
xlab="Chromosome", ylab=expression(-log[10](p-value)),
col=c("gray","darkgray"), panel.extra=NULL, pch=20, cex=0.8,...) {

if (length(chr)==0) stop("chromosome vector is empty")
if (length(pos)==0) stop("position vector is empty")
if (length(pvalue)==0) stop("pvalue vector is empty")

#make sure we have an ordered factor
if(!is.ordered(chr)) {
chr <- ordered(chr)
} else {
chr <- chr[,drop=T]
}

#make sure positions are in kbp
if (any(pos>1e6)) pos<-pos/1e6;

#calculate absolute genomic position
#from relative chromosomal positions
posmin <- tapply(pos,chr, min);
posmax <- tapply(pos,chr, max);
posshift <- head(c(0,cumsum(posmax)),-1);
names(posshift) <- levels(chr)
genpos <- pos + posshift[chr];
getGenPos<-function(cchr, cpos) {
p<-posshift[as.character(cchr)]+cpos
return(p)
}

#parse annotations
grp <- NULL
ann.settings <- list()
label.default<-list(x="peak",y="peak",adj=NULL, pos=3, offset=0.5,
col=NULL, fontface=NULL, fontsize=NULL, show=F)
parse.label<-function(rawval, groupname) {
r<-list(text=groupname)
if(is.logical(rawval)) {
if(!rawval) {r$show <- F}
} else if (is.character(rawval) || is.expression(rawval)) {
if(nchar(rawval)>=1) {
r$text <- rawval
}
} else if (is.list(rawval)) {
r <- modifyList(r, rawval)
}
return(r)
}

if(!is.null(annotate)) {
if (is.list(annotate)) {
grp <- annotate[[1]]
} else {
grp <- annotate
}
if (!is.factor(grp)) {
grp <- factor(grp)
}
} else {
grp <- factor(rep(1, times=length(pvalue)))
}

ann.settings<-vector("list", length(levels(grp)))
ann.settings[[1]]<-list(pch=pch, col=col, cex=cex, fill=col, label=label.default)

if (length(ann.settings)>1) {
lcols<-trellis.par.get("superpose.symbol")$col
lfills<-trellis.par.get("superpose.symbol")$fill
for(i in 2:length(levels(grp))) {
ann.settings[[i]]<-list(pch=pch,
col=lcols[(i-2) %% length(lcols) +1 ],
fill=lfills[(i-2) %% length(lfills) +1 ],
cex=cex, label=label.default);
ann.settings[[i]]$label$show <- T
}
names(ann.settings)<-levels(grp)
}
for(i in 1:length(ann.settings)) {
if (i>1) {ann.settings[[i]] <- modifyList(ann.settings[[i]], ann.default)}
ann.settings[[i]]$label <- modifyList(ann.settings[[i]]$label,
parse.label(ann.settings[[i]]$label, levels(grp)[i]))
}
if(is.list(annotate) && length(annotate)>1) {
user.cols <- 2:length(annotate)
ann.cols <- c()
if(!is.null(names(annotate[-1])) && all(names(annotate[-1])!="")) {
ann.cols<-match(names(annotate)[-1], names(ann.settings))
} else {
ann.cols<-user.cols-1
}
for(i in seq_along(user.cols)) {
if(!is.null(annotate[[user.cols[i]]]$label)) {
annotate[[user.cols[i]]]$label<-parse.label(annotate[[user.cols[i]]]$label,
levels(grp)[ann.cols[i]])
}
ann.settings[[ann.cols[i]]]<-modifyList(ann.settings[[ann.cols[i]]],
annotate[[user.cols[i]]])
}
}
rm(annotate)

#reduce number of points plotted
if(should.thin) {
thinned <- unique(data.frame(
logp=round(-log10(pvalue),thin.logp.places),
pos=round(genpos,thin.pos.places),
chr=chr,
grp=grp)
)
logp <- thinned$logp
genpos <- thinned$pos
chr <- thinned$chr
grp <- thinned$grp
rm(thinned)
} else {
logp <- -log10(pvalue)
}
rm(pos, pvalue)
gc()

#custom axis to print chromosome names
axis.chr <- function(side,...) {
if(side=="bottom") {
panel.axis(side=side, outside=T,
at=((posmax+posmin)/2+posshift),
labels=levels(chr),
ticks=F, rot=0,
check.overlap=F
)
} else if (side=="top" || side=="right") {
panel.axis(side=side, draw.labels=F, ticks=F);
}
else {
axis.default(side=side,...);
}
}

#make sure the y-lim covers the range (plus a bit more to look nice)
prepanel.chr<-function(x,y,...) {
A<-list();
maxy<-ceiling(max(y, ifelse(!is.na(sig.level), -log10(sig.level), 0)))+.5;
A$ylim=c(0,maxy);
A;
}

xyplot(logp~genpos, chr=chr, groups=grp,
axis=axis.chr, ann.settings=ann.settings,
prepanel=prepanel.chr, scales=list(axs="i"),
panel=function(x, y, ..., getgenpos) {
if(!is.na(sig.level)) {
#add significance line (if requested)
panel.abline(h=-log10(sig.level), lty=2);
}
panel.superpose(x, y, ..., getgenpos=getgenpos);
if(!is.null(panel.extra)) {
panel.extra(x,y, getgenpos, ...)
}
},
panel.groups = function(x,y,..., subscripts, group.number) {
A<-list(...)
#allow for different annotation settings
gs <- ann.settings[[group.number]]
A$col.symbol <- gs$col[(as.numeric(chr[subscripts])-1) %% length(gs$col) + 1]
A$cex <- gs$cex[(as.numeric(chr[subscripts])-1) %% length(gs$cex) + 1]
A$pch <- gs$pch[(as.numeric(chr[subscripts])-1) %% length(gs$pch) + 1]
A$fill <- gs$fill[(as.numeric(chr[subscripts])-1) %% length(gs$fill) + 1]
A$x <- x
A$y <- y
do.call("panel.xyplot", A)
#draw labels (if requested)
if(gs$label$show) {
gt<-gs$label
names(gt)[which(names(gt)=="text")]<-"labels"
gt$show<-NULL
if(is.character(gt$x) | is.character(gt$y)) {
peak = which.max(y)
center = mean(range(x))
if (is.character(gt$x)) {
if(gt$x=="peak") {gt$x<-x[peak]}
if(gt$x=="center") {gt$x<-center}
}
if (is.character(gt$y)) {
if(gt$y=="peak") {gt$y<-y[peak]}
}
}
if(is.list(gt$x)) {
gt$x<-A$getgenpos(gt$x[[1]],gt$x[[2]])
}
do.call("panel.text", gt)
}
},
xlab=xlab, ylab=ylab,
panel.extra=panel.extra, getgenpos=getGenPos, ...
);
}
</syntaxhighlight>

= Sample Usage =

First, let us create some fake data to test it with. We will use an imaginary animal that has 10 autosomes. The random sample data is created with
<syntaxhighlight lang="rsplus">
#FAKE SAMPLE DATA
createSampleGwasData<-function(chr.count=10, include.X=F) {
chr<-c(); pos<-c()
for(i in 1:chr.count) {
chr <- c(chr,rep(i, 1000))
pos <- c(pos,ceiling(runif(1000)*(chr.count-i+1)*25*1e3))
}
if(include.X) {
chr <- c(chr,rep("X", 1000))
pos <- c(pos,ceiling(runif(1000)*5*25*1e3))
}
pvalue <- runif(length(pos))
return(data.frame(chr, pos,pvalue))
}
dd<-createSampleGwasData()
dd$pvalue[3000] <- 1e-7 #include a significant result
</syntaxhighlight>
Using this sample data (chr, pos, and pvalue) we can make the manhattan plot by calling
<syntaxhighlight lang="rsplus">
manhattan.plot(dd$chr, dd$pos, dd$pvalue)
</syntaxhighlight>
[[File:Sample-manhattan-plot-A.png|center]]
That's it! But if that were all the function could do, it would be much shorter. There are many other parts of the plot that you can customize.

The <tt>sig.level=</tt> parameter allows you to draw a dotted line indicating the genome-wide significance threshold for your data. If you leave out that parameter, no line will be drawn. A common option would be to incluide <tt>sig.level=5e-8</tt> for a genome-wide association study.

Notice how the colors alternate between chromosomes. You can change these colors with the <tt>col=</tt> parameter. You need to pass in a vector of R colors. For example, if you want a more festive plot, try <tt>col=c("orange","blue","purple")</tt>.

By default the function attempts to minimize the number of points drawn by rounding the -log10 p-value and the position and then only plotting the unique combinations. You can disable this default behavior by specifying <tt>should.thin=F</tt> when you call the function. You can also contol the level of thinning by setting <tt>thin.pos.places=</tt> and <tt>thin.logp.places=</tt>. These values control the number of decimal places that each value is rounded to. Both of these values default to 2 decimal places.

Furthermore, you are able to pass any additional lattice graphical parameters that you want as well such as <tt>main=</tt>, <tt>ylim=</tt>, or <tt>par.settings=</tt>. They will be passed though to the embded <tt>xyplot()</tt> command.

If your list of chromosomes is non-numeric, you need to bit of extra work. You must make sure the value you pass in for chromosome is an ordered factor, such that chromosomes 1-22 come first, followed by X,Y, and/or MT. An example might be
<syntaxhighlight lang="rsplus">
#generate sample data (with X) using above function
dd<-createSampleGwasData(include.X=T)
#now plot chromosomes in correct order
manhattan.plot(factor(dd$chr, levels=c(1:10, "X")), dd$pos, dd$pvalue)
</syntaxhighlight>
Noteice that the labels on the x-axis in the plot come from the levels of the chr factor. Thus if you have raw data that has chromosomes 1-25 where stands for 23=X, 24,Y, and 25=MT, you can create the appropriate ordered vector with the proper names using <tt>factor(chr, levels=1:25, labels=c(1:22, "X","Y","MT"))</tt> where <tt>chr</tt> is your vector of values 1-25.

In tests, running R to read in GWAS results (2.5 million SNPs) and create a manhattan plot using this function took about 7-10 minutes. The only real concern is how much memory R uses when you read in the data. It is important to only read in the data that you need for the plot to minimize memory; so if your results file contains other columns, you may wish to ignote them using a NULL in <tt>colClasses=</tt> (see the documentation for <tt>read.table</tt> for more information).

= Annotating Regions or SNPs =

The <tt>mantahhan.plot</tt> function provies many options for annotating differnt parts of your plot. For example you may wish to highlight certain gene regions or point out certain SNPs. You can do this with the <tt>annotate=</tt> parameter. In the simplest case, you can pass in a factor (with the same length as the pvalue vector) which assigns each point to a group. The function assumes the '''first level''' of that factor correspondes to the background SNPs, and every other factor will be drawn in a different color, and labeled with their lable name. (Note: you can check the levels and labels of a factor variable in R with <tt>levels()</tt> and <tt>labels()</tt>). Here is an example of coloring in three gene regions
<syntaxhighlight lang="rsplus">
#create sample data
dd<-createSampleGwasData()
#make annotation factor
ann<-rep(1, length(dd$pvalue))
ann[with(dd, chr==1 & pos>=90e3 & pos<110e3)]<-2
ann[with(dd, chr==4 & pos>=50e3 & pos<80e3)]<-3
ann[with(dd, chr==6 & pos>=30e3 & pos<50e3)]<-4
ann<-factor(ann, levels=1:4, labels=c("","GENE1","GENE2", "GENE3"))
#draw plot with annotation
manhattan.plot(dd$chr, dd$pos, dd$pvalue, annotate=ann)
</syntaxhighlight>

[[File:Sample-manhattan-plot-B.png|center]]

Notice that the new regions (those other than the first level of the factor) are colored accooring to the default trellis superpose.symbol colors. However, you have a lot of control over how these special regions are drawn and labeled. Rather than passing in a factor to <tt>annotate=</tt>, you can pass in a list containing your factor as well as additional lists to customize the drawing of each of the different regions. Each region has the following settings:
<syntaxhighlight lang="rsplus">
list(col= pch=, cex=, fill=,
label=list(text=, x=,y=,adj=, pos=, offset=,
col=, fontface=, fontsize=, show=)
)
</syntaxhighlight>
All of these parmeters have there usual meaning for plotting functions (see <tt>?xyplot</tt> or <tt>?ltext</tt> for more information on these parameters). So when building a list for your annotation,
your factor must come first, then you can also pass in lists corresponding to each level of your
factor (including the background SNPs) that contain the values you want to override. Alternatively, you can pass in lists with names corresponding to the labels of your factor. Let's illustrate with some examples. Say we wanted to change the plotting character of GENE2 to be trangles, we can either do
<syntaxhighlight lang="rsplus">
manhattan.plot(dd$chr, dd$pos, dd$pvalue,
annotate=list(ann, list(), list(), list(pch=17), list()))
#or
manhattan.plot(dd$chr, dd$pos, dd$pvalue,
annotate=list(ann, "GENE2"=list(pch=17)))
</syntaxhighlight>
Notice in the first example, there are four <tt>list()</tt> items in the annotate list, each corresponse to the four levels of the <tt>ann</tt> factor. In the second example, we specifically set properties for just the "GENE2" level.

So we can also customize how the labels are displayed. You can change the text with any of the following parameters
<syntaxhighlight lang="rsplus">
annotate=list(ann, "GENE2"=list(label="Other Name"))
#or
annotate=list(ann, "GENE2"=list(label=list(text="Other Name")))
</syntaxhighlight>
or you can turn off labeling with either of the following
<syntaxhighlight lang="rsplus">
annotate=list(ann, "GENE2"=list(label=F))
#or
annotate=list(ann, "GENE2"=list(label=list(show=F)))
</syntaxhighlight>
Further more, you can specify an <tt>x=</tt> and <tt>y=</tt> value ```in a label list``` to control where the label is placed. For <tt>x=</tt>, you can speicfy an chromosome and position in the format <tt>x=list(chr,pos)</tt> or you can specify "center" to find the x-value that's in the exact middle of the region, or "peak" to get the x-value at the SNP with the most significant p-value. The default is "peak". For the <tt>y=</tt> parameter, you can specify a -log10 pvalue (not just the regualr pvalue) or "peak" to get the -log10 pvalue for the best SNP in the region. All of these are valid locations:
<syntaxhighlight lang="rsplus">
annotate=list( ann, "GENE2"=list(label=list(x="center")) )
annotate=list( ann, "GENE2"=list(label=list(x="peak", y=7)) )
annotate=list( ann, "GENE2"=list(label=list(x=list(3,5e3))) )
</syntaxhighlight>
The default position for the label, in reference to the (x,y) point, is <tt>pos=3</tt> which indicates the label is to be centered above the point. If you want this to be even a bit higher than than point, you can also increase the offset with <tt>offset=2</tt> or something similar. For more info on these options, see <tt>?text</tt>. Here is what the list might locate if we wanted to change the color of a gene region and mvoe the label somewhat:
<syntaxhighlight lang="rsplus">
annotate=list( ann, "GENE2"=list(col="green", label=list(x="center", offset=2)) )
</syntaxhighlight>

The defaults for the first level are taken from the <tt>pch=</tt> and <tt>col=</tt> parameters you pass to the function. If you want to change the defaults for the additional gene regions, you can set the <tt>ann.default=</tt> parameter. These are the settings applied to the regions before any of the settings you might specify in the <tt>annotate=</tt> parameter. For example, if you wanted to make all the highlighed regions use triangles, rather than circles, you could run
<syntaxhighlight lang="rsplus">
manhattan.plot(dd$chr, dd$pos, dd$pvalue, annotate=ann, ann.default=list(pch=17))
</syntaxhighlight>

= Annotation Helper Function =

While you have a lot of options to put together an annotation list, there are some types of annotations that are quite standard. To help with these cases, you may use this <tt>annotateSNPRegions()</tt> function to easily create annotations around SNPs.

<syntaxhighlight lang="rsplus">
annotateSNPRegions<-function(snps, chr, pos, pvalue, snplist,
kbaway=0, maxpvalue=1, labels=c(), col=c(), pch=c()) {

stopifnot(all(length(snps)==length(chr), length(chr)==length(pos),
length(pos)==length(pvalue)))
if (length(snplist)==0) stop("snplist vector is empty")

if(any(pos>1e6)) kbaway<-kbaway*1000

ann<-rep(0, length(snps))
for(i in seq_along(snplist)) {
si<-which(snps==snplist[i])
ci<-chr[si]
pi<-pos[si]
ann[chr==ci & pos >= pi-kbaway & pos <= pi+kbaway & pvalue<=maxpvalue]<-i
}
ann<-list(factor(ann, levels=0:length(snplist), labels=c("", snplist)))
if(length(col)>0 || length(pch)>0 || length(labels)>0) {
for(i in seq_along(snplist)) {
ann[[ snplist[i] ]] = list()
if(length(col)>0) {
ann[[ snplist[i] ]]$col = col[ (i-1) %% length(col)+1 ]
}
if(length(pch)>0) {
ann[[ snplist[i] ]]$pch = pch[ (i-1) %% length(pch)+1 ]
}
if(length(labels)>0) {
ann[[ snplist[i] ]]$label = labels[ (i-1) %% length(labels)+1 ]
}
}
}
return(ann)
}
</syntaxhighlight>

To use this function, you must pass in snp names, chromosome, position, and pvalue. Then you must specify a list of SNPs that you want to use to identify your regions. Then you can optionally specify a color (<tt>col=</tt>) or plotting character (<tt>pch=</tt>) to use for those regions. If you give more than one color (for example), those colors will be cycled over all of the regions. If you supply only one color, the same color will be used for all regions

We can show how to use this function with sample sample data. Here is an example where we identiy 3 gene regions:
<syntaxhighlight lang="rsplus">
#generate fake data
dd<-createSampleGwasData()
dd$snp <- paste("rs", 1:length(dd$pvalue), sep="") #make up SNP names
#create annotations
ann<-annotateSNPRegions(dd$snp, dd$chr,dd$pos,dd$pvalue,
c("rs10", "rs1001", "rs2005"),
labels=c("GENE1","GENE2","GENE3"),
col="blue",
kbaway=100
)
#draw plot
manhattan.plot(dd$chr, dd$pos, dd$pvalue, annotate=ann)
</syntaxhighlight>

= A final example =

Here is the code used to create the image at the top of the page

<syntaxhighlight lang="rsplus">
dd<-read.table("/path/to/data.txt", header=T, as.is=T, sep="\t")
ann<-annotateSNPRegions(dd$MarkerName, dd$CHR, dd$POS, dd$P.value,
snplist=c("rs10994397","rs7296288","rs736408","rs9371601", "rs12576775"),
name=c("ANK3","CHR12","ITIH3","SYNE1","ODZ4"),
col=c("green","red")[c(1,1,1,2,2)], kbaway=50
)
png("mh.png", width=950, height=400)
print(manhattan.plot(dd$CHR, dd$POS, dd$P.value,
annotate=ann, ann.default=list(label=list(offset=2)),
sig.level=5e-8,
key=list(background="white", border=T, padding.text=2,
corner=c(.95, .95), text=list(lab=c("Known","Novel")),
points=list(col=c("red","green"), pch=20))
)
)
dev.off()
</syntaxhighlight>

Note that we use the built-in lattice parameter <tt>key=</tt> to create a color key for our gene regions. We cannot use <tt>auto.key=</tt> or <tt>simpleKey()</tt> since this technique uses groups differently than standard Lattice does in order to provide greater flexibility.

= Creating Image Files =

If you have a lot of data, drawing the image to the screen may take a long time, even with the default thinning. Similarly, if you write the image to a PDF, you will get a larger image size and it may take a long time to open/print since the location of every point is stored in the PDF. If you create a PNG file, the image size will be a lot smaller since it remembers the final pixels of the image rather than all the points it drew to make those pixels. If you make the image large enough, you probably won't notice much loss of detail either. Something like
<syntaxhighlight lang="rsplus">
png("mymanhattan.png", width=950, height=500)
print(manhattan.plot(chr, pos, pvalue, sig.level=5e-8))
dev.off()
</syntaxhighlight>
should work well.

= Alternatives =

Note that Cristen Willer and Serena Sanna have also made many beautiful plots for GWAS data using R base graphics. They have provided a sample script that does many of the same things as the above function, plus takes care of reading in the data, highlighting known SNPs, and making QQ plots. The sample code can be found at
snowwhite:/home/cristen/scripts/makeplot_manhattan_qq_forsharing.R

[[Category:Code Samples]]

Code Sample: Generating Manhattan Plots in R - Revision history

Mflick at 14:49, 14 January 2014

Goncalo at 19:53, 9 January 2013

Goncalo: Created page with 'A useful way to summarize genome-wide association data is with a Manhattan plot. This type of plot has a point for every SNP or location tested with the position in the genome al…'