```r
# # Here - use a smaller data set, already filtered for min coverage 5x 
# meth_unite_5x <- methylKit::unite(meth_5x, destrand=TRUE, min.per.group=1L)
# save(meth_unite_5x, file = "../results/methylkit/meth_unite_5x") #save object to file
# load(file = "../results/methylkit/meth_unite_5x") #save object to file
# 
# meth_unite_5x_reshaped <- meth_unite_5x %>%
#   melt(id=c("chr", "start", "end", "strand"), value.name = "count") %>%
#   drop_na(count) %>%
#   mutate(stat=gsub('[0-9]+', '', variable)) %>%
#   mutate(sample=gsub('[a-zA-Z]+', '', variable)) %>%
#   dplyr::select(-variable) %>%
#   pivot_wider(names_from = stat, values_from = count) %>%
#   mutate(percMeth = 100*(numCs/coverage)) %>%
#   mutate(sample=as.numeric(sample))
# 
# save(meth_unite_5x_reshaped, file = "meth_unite_5x_reshaped")
# save(meth_unite_5x_reshaped, file = "../results/methylkit/meth_unite_5x_reshaped")
```