Stat Manhattan

stat_manhattan.R

# # MIT License
#
# Copyright (c) Mickaël Canouil
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# library(data.table)
# library(ggplot2)
# library(scales)

fortify.manhattan <- function(data, x, y, group) {
  map_chro <- c(seq(22), "X", "Y", "X", "Y")
  names(map_chro) <- c(seq(24), "X", "Y")
  
  `:=` <- data.table::`:=`
  
  out <- data.table::as.data.table(data)
  data.table::setnames(out, c(x, y, group), c("x_pos", "y_pval", "x_chr"))
  out[, x_chr := as.character(x_chr)]
  out[, x_chr := map_chro[gsub("^chr", "", x_chr, ignore.case = TRUE)]]
  out[, x_chr := factor(x_chr, levels = intersect(c(seq(22), "X", "Y"), x_chr))]
  out[, x_pos := as.double(x_pos)]
  out[order(x_chr, x_pos)]
  out[, x_pos := scales::rescale(x = x_pos, to = c(-0.4, 0.4)), by = "x_chr"]
  out[, x_pos := x_pos + as.integer(x_chr)]
  data.table::setnames(out, c("x_pos", "y_pval", "x_chr"), c("x", "y", "group"))
  out[]
}

StatManhattan <- ggplot2::ggproto("StatManhattan", ggplot2::Stat,
  required_aes = c("x", "y", "group"),
  setup_data = function(data, params) {
    fortify.manhattan(data, "x", "y", "group")
  },
  compute_layer = function(data, scales, params) {
    data
  }
)

library(ggplot2)
#> Warning: package 'ggplot2' was built under R version 4.0.2
library(scales)
library(data.table)
devtools::source_gist("d7cbc59a7fcbb6d231f432801ec7aa19") # pval_trans
#> Sourcing https://gist.githubusercontent.com/mcanouil/d7cbc59a7fcbb6d231f432801ec7aa19/raw/ae004b3be9110ebb92fa623b8d52076f2b223ac0/pval_trans.R
#> SHA-1 hash of file is 59a871b85174513bf8ee12ca1be9ce4d69f73305
devtools::source_gist("c519e2ade2614f5a375906c3485c4eea") # stat_manhattan
#> Sourcing https://gist.githubusercontent.com/mcanouil/c519e2ade2614f5a375906c3485c4eea/raw/839c2e68ae9e985ecd3c01e6ec18d0b47de03439/stat_manhattan.R
#> SHA-1 hash of file is 0b571cd7549da71fcd13542a220cb9f65a6fbbcc

p <- ggplot(
  data = data.table(
    chr = rep(c(seq(22), "X", "Y"), each = 10),
    pos = 1:240,
    pvalue = rep(10^-c(1:24), each = 10)
  )
) +
  aes(x = pos, y = pvalue) +
  geom_point(stat = "manhattan", size = 0.60, na.rm = TRUE) +
  scale_y_continuous(trans = "pval", expand = expansion(mult = c(0, 0.2)), limits = c(1, NA)) +
  scale_colour_viridis_d(name = "Chr")

p + aes(colour = chr)

p + aes(colour = factor(chr, levels = c(seq(22), "X", "Y")))

library(ggplot2)
#> Warning: package 'ggplot2' was built under R version 4.0.2
library(scales)
library(data.table)
devtools::source_gist("d7cbc59a7fcbb6d231f432801ec7aa19") # pval_trans
#> Sourcing https://gist.githubusercontent.com/mcanouil/d7cbc59a7fcbb6d231f432801ec7aa19/raw/ae004b3be9110ebb92fa623b8d52076f2b223ac0/pval_trans.R
#> SHA-1 hash of file is 59a871b85174513bf8ee12ca1be9ce4d69f73305
devtools::source_gist("c519e2ade2614f5a375906c3485c4eea") # stat_manhattan
#> Sourcing https://gist.githubusercontent.com/mcanouil/c519e2ade2614f5a375906c3485c4eea/raw/839c2e68ae9e985ecd3c01e6ec18d0b47de03439/stat_manhattan.R
#> SHA-1 hash of file is 0b571cd7549da71fcd13542a220cb9f65a6fbbcc

ggplot(
  data = fortify.manhattan(
    data = data.table(
      chr = rep(c(seq(22), "X", "Y"), each = 10),
      pos = 1:240,
      pvalue = rep(10^-c(1:24), each = 10)
    ), 
    x = "pos", 
    y = "pvalue", 
    group = "chr"
  )
) +
  aes(x = x, y = y, colour = group) +
  geom_point(size = 0.60, na.rm = TRUE) +
  scale_y_continuous(trans = "pval", expand = expansion(mult = c(0, 0.2)), limits = c(1, NA))