tidy_matching.R

library(plyranges)
library(readr)

bindata <- read_tsv("bindata.40000.hg19.tsv.gz")

fire <- read_bed("fire-adult-hg19.txt")

# not necessary, but nice to have
si <- Seqinfo(genome="hg19")
si <- keepStandardChromosomes(si)
seqlevels(fire) <- seqlevels(si)
seqinfo(fire) <- si

x <- bindata |>
  select(seqnames=chr, end, fracN, fracGC, zooBasesFdr05) |>
  mutate(width = 40e3) |>
  as_granges()

seqlevels(x) <- seqlevels(si)
seqinfo(x) <- si

# check expected overlap
length(fire)
x |>
  join_overlap_inner(fire) |>
  length()

# filter NAs
x <- x |>
  filter(!is.na(fracGC) & !is.na(zooBasesFdr05))

# don't ask why but this pipe has to be %>%
x <- x %>%
  mutate(fire = count_overlaps(., fire))

table(x$fire)

library(ggplot2)
library(tibble)
library(tidyr)
x |>
  as_tibble() |>
  pivot_longer(cols=c("fracN","fracGC","zooBasesFdr05")) |>
  mutate(fire = factor(fire)) |>
  ggplot(aes(value, fill=fire)) +
  geom_histogram(aes(y = after_stat(density)), position="dodge") +
  facet_wrap(~name, scales="free")

library(nullranges)
set.seed(123)
m <- matchRanges(focal = x %>% filter(fire == 1),
                 pool = x %>% filter(fire == 0),
                 covar = ~fracGC)

plotCovariate(m)

# bad variable naming Mike!
y <- bind_ranges(
  focal = focal(m),
  matched = matched(m),
  .id="type"
)

dat <- y |>
  as_tibble() |>
  mutate(fire = factor(fire))

dat |>
  ggplot(aes(zooBasesFdr05, fill=fire)) +
  geom_histogram(position="dodge")

dat |>
  filter(zooBasesFdr05 < 5000) |>
  ggplot(aes(fire, zooBasesFdr05)) +
  geom_violin() +
  stat_summary(fun = "mean",
               geom = "point",
               colour = "red")

zb <- split(dat$zooBasesFdr05, dat$fire)
t.test(zb[["0"]], zb[["1"]])

library(lsr)
cohensD(zb[["0"]], zb[["1"]])