0

I am creating a Manhattan plot in R, where my aim is to highlight SNPs based on the trait they relate to.

Here is the tail of my main dataframe (gwas_plot):

              SNP CHR        BP         P 
58008  rs77855001   1 249213034 8.665e-01
58009  rs74322942   1 249218520 7.437e-01
58010 rs114152373   1 249222547 6.504e-02
58011  rs11205304   1 149906423 4.000e-23
58012  rs11205305   1 149906433 2.000e-16
58013   rs1061956   1 149914327 4.000e-06

I have another dataset (gwascat_topsnp) which contains a subset of the SNPs in gwas_plot, with a description included. It looks like this:

         SNP CHR        BP     P            Description
1 rs11205301   1 149906412 4e-23                 Height
2 rs11205301   1 149906412 2e-16                 Height
3  rs1061954   1 149914353 4e-06 Obesity-related Traits

I altered the source code for the manhattan plot (original source code found here: https://github.com/stephenturner/qqman/tree/master/R), at a particular section which allows you to highlight specific SNPs. I altered it so that the pch would be unique for each "Description" category (in the above data frame).

This is the section I added:

  keyword2=gwascat_topsnp$Description
# Highlight2: highlight snps from a second character vector
  if (!is.null(highlight2)) {
    if (any(!(highlight2 %in% d$SNP))) warning("You're trying to highlight SNPs that don't exist in your results.")
    d.highlight2=d[which(d$SNP %in% highlight2), ]
    with(d.highlight2, points(pos, logp, col="red1", pch=1:24[as.factor(keyword2)], ...)) 
  }

So here is the code that I use to plot the Manhattan plot:

par(xpd=T, mar=par()$mar+c(0,0,3,0))
manhattan_KB(gwas_plot,col=c("grey"),highlight=gwas_main[rownum,1],highlight2=gwascat_topsnp$SNP)
    par(xpd=TRUE)
    legend(0,40,legend=levels(as.factor(keyword2)),
           col=c("red1"),pch=1:24[as.factor(keyword2)],
           bty="n",cex=0.8,ncol=1,horiz=F)
    par(mar=c(5,4,4,2) + 0.1)

However, in the output seen here, both SNPs with the description labelled as "Height" should have the same pch, however they do not: One of them has a circle and one of them has a triangle, and the other trait remians unlabelled in the legend.

Could someone please help me with where I am going wrong? If you need any more information I will be happy to provide it.

Thank you!

IcedCoffee
  • 375
  • 2
  • 14
  • Can you make your example [reproducible](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example)? – Roman Luštrik Feb 09 '16 at 19:00
  • Do you get any warnings? Would appear that the code `1:24[as.factor(keyword2)]` is not valid and produces a warning. When I test this with the following I am able to reproduce the issue partially: `plot(c(1,2,3), c(1,2,3), pch=1:24[as.factor(c("Height","Height","Obesity"))])` – Vince Feb 22 '16 at 21:08

0 Answers0