Custom manhattan plot multi x-axis

279 Views Asked by At

I have the following data set gwas_data

Running head -n 23 gwas_data gives me the following table.

gwas_data <- 
  
  data.frame(
  stringsAsFactors = FALSE,
                 udi = c("A","B","C","D","E",
                         "F","G","H","I","J","K","A","B","C","D","E",
                         "F","G","H","I","J","K"),
                 snp = c("rs71628639_A",
                         "rs71628639_A","rs71628639_A","rs71628639_A","rs71628639_A",
                         "rs71628639_A","rs71628639_A","rs71628639_A",
                         "rs71628639_A","rs71628639_A","rs71628639_A","rs12726330_A",
                         "rs12726330_A","rs12726330_A","rs12726330_A",
                         "rs12726330_A","rs12726330_A","rs12726330_A","rs12726330_A",
                         "rs12726330_A","rs12726330_A","rs12726330_A"),
                 chr = c(1L,1L,1L,1L,1L,1L,1L,
                         1L,1L,1L,1L,1L,1L,1L,1L,1L,1L,1L,1L,1L,1L,
                         1L),
                  bp = c(154988255L,154988255L,
                         154988255L,154988255L,154988255L,154988255L,154988255L,
                         154988255L,154988255L,154988255L,154988255L,
                         155108167L,155108167L,155108167L,155108167L,155108167L,
                         155108167L,155108167L,155108167L,155108167L,
                         155108167L,155108167L),
                   p = c(0.580621191,0.356577427,
                         0.494774059,0.984005886,0.492034614,0.581479389,
                         0.24820214,0.202720896,0.295462221,0.845848783,
                         0.954714162,0.343101621,0.740942238,0.929127071,0.717965027,
                         0.335111376,0.857154424,0.480087195,0.980307843,
                         0.521114038,0.583150471,0.925783695),
                beta = c(0.000852277,0.003943912,
                         0.001091986,-3.18e-05,0.000564413,0.000120028,
                         0.026156467,0.000303135,0.069146449,-2.96e-07,-2.11e-05,
                         0.001274261,-0.001232397,0.000123948,-0.000498507,
                         -0.000689988,-3.41e-50,-0.013934416,5.12e-06,
                         -0.03696031,-7.28e-07,-3.01e-05),
              bp_cum = c(1.154988255,1.154988255,
                         1.154988255,1.154988255,1.154988255,1.154988255,
                         1.154988255,1.154988255,1.154988255,1.154988255,
                         1.154988255,1.155108167,1.155108167,1.155108167,
                         1.155108167,1.155108167,1.155108167,1.155108167,1.155108167,
                         1.155108167,1.155108167,1.155108167)
  )

I would like to make a manhattan plot, the X-axis should have chromosomal numbers from 1:22, I want each entry to be on the x-axis according to the BP position. The id should act as colour and the y-axis would be -log10(p).

I have rewritten the r command as follows, but my graph doesn't look correct.

library(plyr)
library(dplyr)
library(purrr)
library(tidyverse)
library(ggtext)
library(stringr)
gwas_data <- read.table("gwas_data", header=T)
sig <- 5e-8
manhplot <- ggplot(gwas_data, aes(x = bp_cum, y = -log10(p), color = udi)) +
   geom_hline(yintercept = -log10(sig), color = "grey40", linetype = "dashed") + 
   geom_point(aes(color=as.factor(udi)), alpha=0.8, size=2) +
   scale_x_continuous(label = axis_set$chr, breaks = axis_set$center) +
   scale_y_continuous(expand = c(0,0), limits = c(0, ylim)) +
   #scale_color_manual(values = rep(c("#276FBF", "#183059"), (length(axis_set$chr)))) +
   scale_size_continuous(range = c(0.5,3)) +
   theme_minimal()
print(manhplot)

I would also like to add the name of the ID and SNP if they are above the significant threshold.

My axis_set looks as follows with test data which goes from chromosome 1:4

chr center
1   179641307
2   354697451
3   553030055
4   558565909

My final graph looks as follows: GWAS

0

There are 0 best solutions below