##### Packages
library(ggrepel)
library(ggplot2)
library(dplyr)

##### Plotting variables
significant <- 5e-8 # significant threshold line
suggestive <- 1e-6 # suggestive threshold line
highlight_snps <- c("rs3057","rs3056","rs3060") # snps to highlight in a different colour to other SNPs
annotate_snps <- c("rs3057","rs3056","rs3060") # snps to annotate with a label of their rsID

##### Manhattan plot function
gg.manhattan <- function(df, CHR, P, BP, threshold, annotate, hlight, col, ylims, title){
  # format df
  df.tmp <- df %>% 
    
    # Compute chromosome size
    group_by(CHR) %>% 
    summarise(chr_len=max(BP)) %>% 
    
    # Calculate cumulative position of each chromosome
    mutate(tot=cumsum(chr_len)-chr_len) %>%
    select(-chr_len) %>%
    
    # Add this info to the initial dataset
    left_join(df, ., by=c("CHR"="CHR")) %>%
    
    # Add a cumulative position of each SNP
    arrange(CHR, BP) %>%
    mutate( BPcum=BP+tot) %>%
    
    # Add highlight and annotation information
    mutate( is_highlight=ifelse(SNP %in% hlight, "yes", "no")) %>%
    mutate( is_annotate=ifelse(P < threshold, "yes", "no")) %>%
    mutate( is_annotate=ifelse(SNP %in% annotate, "yes", "no"))
  
  
  # get chromosome center positions for x-axis
  axisdf <- df.tmp %>% group_by(CHR) %>% summarize(center=( max(BPcum) + min(BPcum) ) / 2 )
  
  ggplot(df.tmp, aes(x=BPcum, y=-log10(P))) +
    # Show all points
    geom_point(aes(color=as.factor(CHR)), alpha=0.8, size=0.5) +
    scale_color_manual(values = rep(col, 22 )) +
    
    # custom X axis:
    scale_x_continuous( label = axisdf$CHR, breaks= axisdf$center ) +
    scale_y_continuous(expand = c(0, 0), limits = ylims) + # expand=c(0,0) removes space between plot area and x axis 
    
    # add plot and axis titles
    ggtitle(paste0(title)) +
    labs(x = "Chromosome") +
    
    # add genome-wide sig and sugg lines
    geom_hline(yintercept = -log10(significant), linetype = "solid", col = "black") +
    geom_hline(yintercept = -log10(suggestive), linetype = "dashed", col = "black") + # comment out if you dont want suggetsive line
    
    # Add highlighted points
    geom_point(data=subset(df.tmp, is_highlight=="yes"), color="black", size=2) +
    
    # Add label using ggrepel to avoid overlapping
     geom_text_repel(data=df.tmp[df.tmp$is_annotate=="yes",],
                     aes(label=SNP),
                     min.segment.length = 0.01,
                     force = 1,
                     nudge_x = 0.5,
                     nudge_y = 3,
                     colour = "black"
                     ) +
   
    
    # Customise the theme:
    theme_bw(base_size = 11) +
    theme( 
      plot.title = element_text(hjust = 0.5),
      legend.position="none",
      panel.border = element_blank(),
      panel.grid.major.x = element_blank(),
      panel.grid.minor.x = element_blank(),
      panel.grid.major.y = element_blank(),
      panel.grid.minor.y = element_blank() 
    )
}

# Plot
gg.manhattan(df = gwas_data, 
             annotate = annotate_snps, # provide a list of SNPs to annotate with labels
             threshold = NA, # what threshold will you highlight SNPs from
             hlight = highlight_snps, # provide a list of SNPs to highlight in a different colour
             col = discrete_wes_pal, # provide colours or a vector of colours
             ylims = c(0,10), # provide minimum and maximum of the Y axis
             title = "") 

LS0tCnRpdGxlOiAiTWFuaGF0dGFuIHBsb3QiCm91dHB1dDoKICBodG1sX2RvY3VtZW50OgogICAgdG9jOiBmYWxzZQogICAgbnVtYmVyX3NlY3Rpb25zOiBmYWxzZQogICAgZGZfcHJpbnQ6IHBhZ2VkCiAgICBjb2RlX2Rvd25sb2FkOiB0cnVlCiAgICBoaWdobGlnaHRlcjogbnVsbAotLS0KCmBgYHtyIHNldHVwLCBldmFsPVRSVUUsIGluY2x1ZGU9RkFMU0UsIGVjaG89RkFMU0UsIGVycm9yPUZBTFNFLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFLCBjYWNoZT1UUlVFLCBmaWcuYWxpZ249J2NlbnRlcicsIGNvbW1lbnQ9IiJ9CmNob29zZUNSQU5taXJyb3IoZ3JhcGhpY3M9RkFMU0UsIGluZD0xMzMpCnNvdXJjZSgic291cmNlL3Bsb3RzLW92ZXJ2aWV3LlIiKQpgYGAKX19fCgpQYWdlIGxhc3QgdXBkYXRlZDogYHIgZm9ybWF0KFN5cy50aW1lKCksICclZCAlQiwgJVknKWAKClwgIAoKIyBFeGFtcGxlcyB7I21hbmhhdHRhbl9wbG90IC50YWJzZXR9CgpGb3IgdGhlc2UgcGxvdHMgd2Ugd2lsbCB1c2UgdGhlIFtHV0FTIGRhdGFdKCNnd2FzX2RhdGEpLgoKYGBge3IgZXZhbD1UUlVFLCBpbmNsdWRlPVRSVUUsIGVjaG89RkFMU0UsIGVycm9yPUZBTFNFLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFLCBjYWNoZT1UUlVFLCBmaWcuYWxpZ249J2NlbnRlcicsIGNvbW1lbnQ9IiJ9CmhlYWQoZ3dhc19kYXRhKQpgYGAKClwgIAoKIyMgZ2cubWFuaGF0dGFuKCkgeyNnZy5tYW5oYXR0YW59ClRoaXMgbWFuaGF0dGFuIHBsb3Qgc2NyaXB0IGlzIHRha2VuIGZyb20gW3BjZ29kZGFyZCBnaXRodWJdKGh0dHBzOi8vZ2l0aHViLmNvbS9wY2dvZGRhcmQvQnVyY2hhcmRsYWJfVHV0b3JpYWxzL3dpa2kvR0dwbG90Mi1NYW5oYXR0YW4tUGxvdC1GdW5jdGlvbikgYmVjYXVzZSBpdCBpcyB3YXkgYmV0dGVyIHRoYW4gd2hhdCBJIGNhbWUgdXAgd2l0aC4KYGBge3IgbWFuaGF0dGFuX3Bsb3QsIGV2YWw9VFJVRSwgaW5jbHVkZT1UUlVFLCBlY2hvPVRSVUUsIGVycm9yPUZBTFNFLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFLCBjYWNoZT1UUlVFLCBmaWcuYWxpZ249J2NlbnRlcicsIGNvbW1lbnQ9IiJ9CiMjIyMjIFBhY2thZ2VzCmxpYnJhcnkoZ2dyZXBlbCkKbGlicmFyeShnZ3Bsb3QyKQpsaWJyYXJ5KGRwbHlyKQoKIyMjIyMgUGxvdHRpbmcgdmFyaWFibGVzCnNpZ25pZmljYW50IDwtIDVlLTggIyBzaWduaWZpY2FudCB0aHJlc2hvbGQgbGluZQpzdWdnZXN0aXZlIDwtIDFlLTYgIyBzdWdnZXN0aXZlIHRocmVzaG9sZCBsaW5lCmhpZ2hsaWdodF9zbnBzIDwtIGMoInJzMzA1NyIsInJzMzA1NiIsInJzMzA2MCIpICMgc25wcyB0byBoaWdobGlnaHQgaW4gYSBkaWZmZXJlbnQgY29sb3VyIHRvIG90aGVyIFNOUHMKYW5ub3RhdGVfc25wcyA8LSBjKCJyczMwNTciLCJyczMwNTYiLCJyczMwNjAiKSAjIHNucHMgdG8gYW5ub3RhdGUgd2l0aCBhIGxhYmVsIG9mIHRoZWlyIHJzSUQKCiMjIyMjIE1hbmhhdHRhbiBwbG90IGZ1bmN0aW9uCmdnLm1hbmhhdHRhbiA8LSBmdW5jdGlvbihkZiwgQ0hSLCBQLCBCUCwgdGhyZXNob2xkLCBhbm5vdGF0ZSwgaGxpZ2h0LCBjb2wsIHlsaW1zLCB0aXRsZSl7CiAgIyBmb3JtYXQgZGYKICBkZi50bXAgPC0gZGYgJT4lIAogICAgCiAgICAjIENvbXB1dGUgY2hyb21vc29tZSBzaXplCiAgICBncm91cF9ieShDSFIpICU+JSAKICAgIHN1bW1hcmlzZShjaHJfbGVuPW1heChCUCkpICU+JSAKICAgIAogICAgIyBDYWxjdWxhdGUgY3VtdWxhdGl2ZSBwb3NpdGlvbiBvZiBlYWNoIGNocm9tb3NvbWUKICAgIG11dGF0ZSh0b3Q9Y3Vtc3VtKGNocl9sZW4pLWNocl9sZW4pICU+JQogICAgc2VsZWN0KC1jaHJfbGVuKSAlPiUKICAgIAogICAgIyBBZGQgdGhpcyBpbmZvIHRvIHRoZSBpbml0aWFsIGRhdGFzZXQKICAgIGxlZnRfam9pbihkZiwgLiwgYnk9YygiQ0hSIj0iQ0hSIikpICU+JQogICAgCiAgICAjIEFkZCBhIGN1bXVsYXRpdmUgcG9zaXRpb24gb2YgZWFjaCBTTlAKICAgIGFycmFuZ2UoQ0hSLCBCUCkgJT4lCiAgICBtdXRhdGUoIEJQY3VtPUJQK3RvdCkgJT4lCiAgICAKICAgICMgQWRkIGhpZ2hsaWdodCBhbmQgYW5ub3RhdGlvbiBpbmZvcm1hdGlvbgogICAgbXV0YXRlKCBpc19oaWdobGlnaHQ9aWZlbHNlKFNOUCAlaW4lIGhsaWdodCwgInllcyIsICJubyIpKSAlPiUKICAgIG11dGF0ZSggaXNfYW5ub3RhdGU9aWZlbHNlKFAgPCB0aHJlc2hvbGQsICJ5ZXMiLCAibm8iKSkgJT4lCiAgICBtdXRhdGUoIGlzX2Fubm90YXRlPWlmZWxzZShTTlAgJWluJSBhbm5vdGF0ZSwgInllcyIsICJubyIpKQogIAogIAogICMgZ2V0IGNocm9tb3NvbWUgY2VudGVyIHBvc2l0aW9ucyBmb3IgeC1heGlzCiAgYXhpc2RmIDwtIGRmLnRtcCAlPiUgZ3JvdXBfYnkoQ0hSKSAlPiUgc3VtbWFyaXplKGNlbnRlcj0oIG1heChCUGN1bSkgKyBtaW4oQlBjdW0pICkgLyAyICkKICAKICBnZ3Bsb3QoZGYudG1wLCBhZXMoeD1CUGN1bSwgeT0tbG9nMTAoUCkpKSArCiAgICAjIFNob3cgYWxsIHBvaW50cwogICAgZ2VvbV9wb2ludChhZXMoY29sb3I9YXMuZmFjdG9yKENIUikpLCBhbHBoYT0wLjgsIHNpemU9MC41KSArCiAgICBzY2FsZV9jb2xvcl9tYW51YWwodmFsdWVzID0gcmVwKGNvbCwgMjIgKSkgKwogICAgCiAgICAjIGN1c3RvbSBYIGF4aXM6CiAgICBzY2FsZV94X2NvbnRpbnVvdXMoIGxhYmVsID0gYXhpc2RmJENIUiwgYnJlYWtzPSBheGlzZGYkY2VudGVyICkgKwogICAgc2NhbGVfeV9jb250aW51b3VzKGV4cGFuZCA9IGMoMCwgMCksIGxpbWl0cyA9IHlsaW1zKSArICMgZXhwYW5kPWMoMCwwKSByZW1vdmVzIHNwYWNlIGJldHdlZW4gcGxvdCBhcmVhIGFuZCB4IGF4aXMgCiAgICAKICAgICMgYWRkIHBsb3QgYW5kIGF4aXMgdGl0bGVzCiAgICBnZ3RpdGxlKHBhc3RlMCh0aXRsZSkpICsKICAgIGxhYnMoeCA9ICJDaHJvbW9zb21lIikgKwogICAgCiAgICAjIGFkZCBnZW5vbWUtd2lkZSBzaWcgYW5kIHN1Z2cgbGluZXMKICAgIGdlb21faGxpbmUoeWludGVyY2VwdCA9IC1sb2cxMChzaWduaWZpY2FudCksIGxpbmV0eXBlID0gInNvbGlkIiwgY29sID0gImJsYWNrIikgKwogICAgZ2VvbV9obGluZSh5aW50ZXJjZXB0ID0gLWxvZzEwKHN1Z2dlc3RpdmUpLCBsaW5ldHlwZSA9ICJkYXNoZWQiLCBjb2wgPSAiYmxhY2siKSArICMgY29tbWVudCBvdXQgaWYgeW91IGRvbnQgd2FudCBzdWdnZXRzaXZlIGxpbmUKICAgIAogICAgIyBBZGQgaGlnaGxpZ2h0ZWQgcG9pbnRzCiAgICBnZW9tX3BvaW50KGRhdGE9c3Vic2V0KGRmLnRtcCwgaXNfaGlnaGxpZ2h0PT0ieWVzIiksIGNvbG9yPSJibGFjayIsIHNpemU9MikgKwogICAgCiAgICAjIEFkZCBsYWJlbCB1c2luZyBnZ3JlcGVsIHRvIGF2b2lkIG92ZXJsYXBwaW5nCiAgICAgZ2VvbV90ZXh0X3JlcGVsKGRhdGE9ZGYudG1wW2RmLnRtcCRpc19hbm5vdGF0ZT09InllcyIsXSwKICAgICAgICAgICAgICAgICAgICAgYWVzKGxhYmVsPVNOUCksCiAgICAgICAgICAgICAgICAgICAgIG1pbi5zZWdtZW50Lmxlbmd0aCA9IDAuMDEsCiAgICAgICAgICAgICAgICAgICAgIGZvcmNlID0gMSwKICAgICAgICAgICAgICAgICAgICAgbnVkZ2VfeCA9IDAuNSwKICAgICAgICAgICAgICAgICAgICAgbnVkZ2VfeSA9IDMsCiAgICAgICAgICAgICAgICAgICAgIGNvbG91ciA9ICJibGFjayIKICAgICAgICAgICAgICAgICAgICAgKSArCiAgIAogICAgCiAgICAjIEN1c3RvbWlzZSB0aGUgdGhlbWU6CiAgICB0aGVtZV9idyhiYXNlX3NpemUgPSAxMSkgKwogICAgdGhlbWUoIAogICAgICBwbG90LnRpdGxlID0gZWxlbWVudF90ZXh0KGhqdXN0ID0gMC41KSwKICAgICAgbGVnZW5kLnBvc2l0aW9uPSJub25lIiwKICAgICAgcGFuZWwuYm9yZGVyID0gZWxlbWVudF9ibGFuaygpLAogICAgICBwYW5lbC5ncmlkLm1ham9yLnggPSBlbGVtZW50X2JsYW5rKCksCiAgICAgIHBhbmVsLmdyaWQubWlub3IueCA9IGVsZW1lbnRfYmxhbmsoKSwKICAgICAgcGFuZWwuZ3JpZC5tYWpvci55ID0gZWxlbWVudF9ibGFuaygpLAogICAgICBwYW5lbC5ncmlkLm1pbm9yLnkgPSBlbGVtZW50X2JsYW5rKCkgCiAgICApCn0KCiMgUGxvdApnZy5tYW5oYXR0YW4oZGYgPSBnd2FzX2RhdGEsIAogICAgICAgICAgICAgYW5ub3RhdGUgPSBhbm5vdGF0ZV9zbnBzLCAjIHByb3ZpZGUgYSBsaXN0IG9mIFNOUHMgdG8gYW5ub3RhdGUgd2l0aCBsYWJlbHMKICAgICAgICAgICAgIHRocmVzaG9sZCA9IE5BLCAjIHdoYXQgdGhyZXNob2xkIHdpbGwgeW91IGhpZ2hsaWdodCBTTlBzIGZyb20KICAgICAgICAgICAgIGhsaWdodCA9IGhpZ2hsaWdodF9zbnBzLCAjIHByb3ZpZGUgYSBsaXN0IG9mIFNOUHMgdG8gaGlnaGxpZ2h0IGluIGEgZGlmZmVyZW50IGNvbG91cgogICAgICAgICAgICAgY29sID0gZGlzY3JldGVfd2VzX3BhbCwgIyBwcm92aWRlIGNvbG91cnMgb3IgYSB2ZWN0b3Igb2YgY29sb3VycwogICAgICAgICAgICAgeWxpbXMgPSBjKDAsMTApLCAjIHByb3ZpZGUgbWluaW11bSBhbmQgbWF4aW11bSBvZiB0aGUgWSBheGlzCiAgICAgICAgICAgICB0aXRsZSA9ICIiKSAKYGBgCgoK