sim_manhattan.R

set.seed(100)
# p1 = runif(1000)
# p2 = runif(5)
z = rep(rnorm(10000))
z[sample(1:10000, size = 3, replace = F)] = rnorm(3, 7, sd = 1)
# plot(-log10(p))
zobs = rep(0, 10000)
density = exp(c(-5:-1, 0, -1:-5)*0.5)
for(i in 6 : 9995) {
  zobs[i] = sum(density * z[c(i-(5:1),i,i+(1:5))])
}
plot(-pnorm(abs(z), lower.tail = F, log.p = T))
df = data.frame(x = -pnorm(abs(zobs), lower.tail = F, log.p = T) -log(2), y = 1:10000)
library(ggplot2)
theme_set(theme_bw(base_size=20))
xinter = c(0, 1700, 3000, 4500, 6000, 8000, 10000)
p = ggplot() + geom_point(data = df, aes(x = y, y = x))
p = p + geom_hline(yintercept = -log(5e-8), linetype = 2) + xlab('genomic position') + ylab('-log(p)') + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black")) + geom_vline(xintercept = xinter, color = 'gray') + coord_cartesian(ylim = c(-30, 50))


pos = c()
lgen = c()
ypos = c()
idx = c()
for(i in 1 : 4) {
  for(j in 1 : (length(xinter) - 1)) {
    ypos = c(ypos, i)
    s = xinter[j]
    e = xinter[j + 1]
    pos = c(pos, runif(1, s, e))
    lgen = c(lgen, runif(1, 100, 200))
    idx = c(idx, j)
  }
}
df2 = data.frame(pos = pos, lgen = lgen, ypos = ypos, idx = idx)

df2$type = 'not selected'
df2$type[df2$idx == 2] = 'selected'
df2$type[df2$idx == 3 & df2$pos < 3500] = 'selected'
df2$is_omim = F
df2$is_omim[c(2, 12)] = T
p = p + geom_errorbarh(data = df2, aes(xmin = pos - lgen, xmax = pos + lgen, y = ypos * -3, color = type), height = 1) + geom_point(data = df2[df2$is_omim, ], aes(x = pos, y = ypos * -3), shape = 17, color = 'red')+ theme(legend.position='none') + scale_color_manual(values = c('not selected' = 'gray', 'selected' = 'black'));p

ggsave('toy_manhattan.png', p, height = 3, width = 6.5)