Intro to the qqman package

The qqman package includes functions for creating manhattan plots and q-q plots from GWAS results. The gwasResults data.frame included with the package has simulated results for 16,470 SNPs on 22 chromosomes. Take a look at the data:

str(gwasResults)
'data.frame':   16470 obs. of  5 variables:
 $ SNP   : chr  "rs1" "rs2" "rs3" "rs4" ...
 $ CHR   : int  1 1 1 1 1 1 1 1 1 1 ...
 $ BP    : int  1 2 3 4 5 6 7 8 9 10 ...
 $ P     : num  0.915 0.937 0.286 0.83 0.642 ...
 $ zscore: num  0.107 0.0789 1.0666 0.2141 0.4653 ...
head(gwasResults)
  SNP CHR BP         P    zscore
1 rs1   1  1 0.9148060 0.1069785
2 rs2   1  2 0.9370754 0.0789462
3 rs3   1  3 0.2861395 1.0666287
4 rs4   1  4 0.8304476 0.2141275
5 rs5   1  5 0.6417455 0.4652597
6 rs6   1  6 0.5190959 0.6447396
tail(gwasResults)
          SNP CHR  BP         P    zscore
16465 rs16465  22 530 0.5643702 0.5763624
16466 rs16466  22 531 0.1382863 1.4822028
16467 rs16467  22 532 0.3936999 0.8529268
16468 rs16468  22 533 0.1778749 1.3473271
16469 rs16469  22 534 0.2393020 1.1767332
16470 rs16470  22 535 0.2630441 1.1192251

How many SNPs on each chromosome?

as.data.frame(table(gwasResults$CHR))
   Var1 Freq
1     1 1500
2     2 1191
3     3 1040
4     4  945
5     5  877
6     6  825
7     7  784
8     8  750
9     9  721
10   10  696
11   11  674
12   12  655
13   13  638
14   14  622
15   15  608
16   16  595
17   17  583
18   18  572
19   19  562
20   20  553
21   21  544
22   22  535

Creating manhattan plots

Now, let's make a basic manhattan plot.

manhattan(gwasResults)

plot of chunk unnamed-chunk-4

We can also pass in other graphical parameters. Let's add a title (main=), increase the y-axis limit (ylim=), reduce the point size to 60% (cex=), and reduce the font size of the axis labels to 90% (cex.axis=). While we're at it, let's change the colors (col=), remove the suggestive and genome-wide significance lines, and supply our own labels for the chromosomes:

manhattan(gwasResults, main = "Manhattan Plot", ylim = c(0, 10), cex = 0.6, 
    cex.axis = 0.9, col = c("blue4", "orange3"), suggestiveline = F, genomewideline = F, 
    chrlabs = c(1:20, "P", "Q"))

plot of chunk unnamed-chunk-5

Now, let's look at a single chromosome:

manhattan(subset(gwasResults, CHR == 1))

plot of chunk unnamed-chunk-6

Let's highlight some SNPs of interest on chromosome 3. The 100 SNPs we're highlighting here are in a character vector called snpsOfInterest. You'll get a warning if you try to highlight SNPs that don't exist.

str(snpsOfInterest)
 chr [1:100] "rs3001" "rs3002" "rs3003" "rs3004" "rs3005" ...
manhattan(gwasResults, highlight = snpsOfInterest)

plot of chunk unnamed-chunk-7

We can combine highlighting and limiting to a single chromosome, and use the xlim graphical parameter to zoom in on a region of interest (between position 200-500):

manhattan(subset(gwasResults, CHR == 3), highlight = snpsOfInterest, xlim = c(200, 
    500), main = "Chr 3")

plot of chunk unnamed-chunk-8

We can also annotate SNPs based on their p-value. By default, this only annotates the top SNP per chromosome that exceeds the annotatePval threshold.

manhattan(gwasResults, annotatePval = 0.01)