Lightweight DNase-seq example

if(!file.exists("data/dat_FDR01_hg38.RData")) {
  download.file(url = "https://zenodo.org/record/3838751/files/dat_FDR01_hg38.RData?download=1",
                destfile = "data/dat_FDR01_hg38.RData",
                method = "curl")
}

if (!file.exists("data/DHS_Index_and_Vocabulary_hg38_WM20190703.txt.gz")) {
  download.file(url = "https://zenodo.org/record/3838751/files/DHS_Index_and_Vocabulary_hg38_WM20190703.txt.gz?download=1",
                destfile = "data/DHS_Index_and_Vocabulary_hg38_WM20190703.txt.gz",
                method = "curl")
  R.utils::gunzip("data/DHS_Index_and_Vocabulary_hg38_WM20190703.txt.gz")
}

if (!file.exists("data/DHS_Index_and_Vocabulary_metadata.xlsx")) {
  download.file(url = "https://zenodo.org/record/3838751/files/DHS_Index_and_Vocabulary_metadata.xlsx?download=1",
                destfile = "data/DHS_Index_and_Vocabulary_metadata.xlsx",
                method = "curl")
}

set.seed(217)

# Extract sample metadata for samples to be used
sample_data <-
  read_xlsx("data/biosamples_used.xlsx",
            range = c("A1:M39"))

# Cell populations to be analyzeds
cell_pops_to_keep <-
  c(
    "CD4.DS17881",
    "CD34.DS12274",
    "CD34_T18.DS25969A",
    "CD14.DS17215",
    "K562.DS16924"
  )


# Join with the data from the source paper's supplement
all_sample_metadata <- 
    read_xlsx("data/DHS_Index_and_Vocabulary_metadata.xlsx") %>%
    right_join(
        sample_data,
        by = c("DCC Library ID" = "DCC_library_id", "DCC Biosample ID" = "DCC_biosample_id")
    )

# Get BED info
bed <- readr::read_tsv("data/DHS_Index_and_Vocabulary_hg38_WM20190703.txt")

# Load in and subset the normalized data (this object is called `dat`)
load("data/dat_FDR01_hg38.RData")

# Subset to selected hematopoietic columns, then filter out rows with 0 or 1 DHSs
my_dat <- dat %>%
    as_tibble() %>%
    dplyr::select(all_sample_metadata$`library order`) %>%
    bind_cols(dplyr::select(bed, seqname, start, end)) %>%
    relocate(seqname, start, end, .before = 1) %>%
    # subset to 5 samples as plotted in the CLIMB article
    dplyr::select(seqname, start, end, matches(cell_pops_to_keep)) %>%
    dplyr::filter(rowSums(across(all_of(4:last_col()), ~ .x != 0)) > 1) %>%
    rename("chr" = "seqname") %>%
    mutate(across(4:last_col(), ~ replace(.x, .x == 0, rnorm(sum(.x==0)))))

# clean large object from environment
rm(dat)

# this is to prevent some numerical issues due to extreme outliers
max_quant <- dplyr::select(my_dat, 4:last_col()) %>%
        map_dbl(~ quantile(.x, .999)) %>%
        max()
        
my_dat <- my_dat %>%
    mutate(across(4:last_col(), ~ replace(.x, .x >= max_quant, max_quant))) %>%
    # select only chromosomes 21 and 22
    filter(chr %in% c("chr21", "chr22")) %>%
    group_split(chr)

for(i in seq_along(my_dat)) {
    if(!dir.exists(paste0("data/DNase_", my_dat[[i]]$chr[1]))) {
        dir.create(paste0("data/DNase_", my_dat[[i]]$chr[1]))
    }
    saveRDS(object = my_dat[[i]], file = paste0("data/DNase_", my_dat[[i]]$chr[1], "/dat.rds"))
}

chr <- 21:22
z <- map(chr, ~ readRDS(paste0("data/DNase_chr", .x, "/dat.rds")) %>%
                 mutate(chr = NULL, start = NULL, end = NULL)) %>%
  set_names(paste0("chr", chr))

fits <- map(z, ~ CLIMB::get_pairwise_fits(z = .x, parallel = TRUE, ncores = 4))

if(!dir.exists("output/DNase")) {
  dir.create("output/DNase")
}

if(!dir.exists("output/DNase/pwfits")) {
  dir.create("output/DNase/pwfits")
}

walk(chr, ~ {
  if (!dir.exists(paste0("output/DNase/pwfits/chr", .x))) {
    dir.create(paste0("output/DNase/pwfits/chr", .x))
  }
})

iwalk(fits, ~ saveRDS(.x, paste0("output/DNase/pwfits/", .y, "/pwfits.rds")))

# This finds the dimension of the data directly from the pairwise fits
D <- as.numeric(strsplit(tail(names(fits[[1]]),1), "_")[[1]][2])

# calculates the sample sizes from the pairwise fits
n <- map_dbl(fits, ~ length(.x[[1]]$cluster))

if(!dir.exists("output/DNase/reduced_classes")) {
  dir.create("output/DNase/reduced_classes")
}
walk(chr, ~ {
  if (!dir.exists(paste0("output/DNase/reduced_classes/chr", .x))) {
    dir.create(paste0("output/DNase/reduced_classes/chr", .x))
  }
})

# Get the list of candidate latent classes
reduced_classes <-
  imap(fits, ~ get_reduced_classes(
    .x,
    D,
    paste0("output/DNase/reduced_classes/", .y, "/lgf.txt"),
    split_in_two = FALSE
  ))

# write the output to a text file
iwalk(reduced_classes, ~ {
  readr::write_tsv(
    data.frame(.x),
    file = paste0("output/DNase/reduced_classes/", .y, "/red_class.txt"),
    col_names = FALSE
  )
})

if(!dir.exists("output/DNase/mcmc")) {
    dir.create("output/DNase/mcmc")
}
walk(chr, ~ {
  if (!dir.exists(paste0("output/DNase/mcmc/chr", .x))) {
    dir.create(paste0("output/DNase/mcmc/chr", .x))
  }
})

# Compute the prior weights
prior_weights <-
  pmap(list(fits, reduced_classes, names(fits)), function(.x, .y, .z)
    get_prior_weights(
      .y,
      .x,
      parallel = FALSE,
      delta = 0:10
    )) %>%
  # just keep all classes since the analysis is small
  map(~ tail(.x, 1)[[1]])
iwalk(prior_weights, ~ saveRDS(.x, paste0("output/DNase/mcmc/", .y, "/prior_weights.rds")))

# obtain the hyperparameters
hyp <-
  pmap(list(my_dat, fits, reduced_classes, prior_weights), function(.w, .x, .y, .z)
    get_hyperparameters(
      as.data.frame(dplyr::select(.w, 4:last_col())),
      .x,
      as.data.frame(.y),
      as.vector(.z)
    )) %>%
  set_names(names(fits))

iwalk(hyp, ~ saveRDS(.x, file = paste0("output/DNase/mcmc/", .y, "/hyperparameters.rds")))

results <-
  pmap(list(my_dat, hyp, reduced_classes), function(.x, .y, .z) run_mcmc(
    dplyr::select(.x, 4:last_col()),
    hyp =  .y,
    nstep = 2000,
    retained_classes = .z
  )) %>%
  set_names(names(fits))

chains <- map(results, extract_chains)

iwalk(chains, ~ saveRDS(.x, file = paste0("output/DNase/mcmc/", .y, "/chain.rds")))

burnin <- 1:500
merged <-
  suppressMessages(merge_classes(
    n_groups = 12,
    # number of classes used in the CLIMB article's analysis
    chain = chains,
    burnin = burnin,
    multichain = TRUE
  ))

col_distmat <- compute_distances_between_conditions(chains, burnin, multichain = TRUE)
row_distmat <- compute_distances_between_clusters(chains, burnin, multichain = TRUE)

colnames(col_distmat) <- names(my_dat[[1]])[-(1:3)]

hc_row <- hclust(as.dist(row_distmat), method = "complete")
hc_col <- hclust(as.dist(col_distmat), method = "complete")

# Get a row reordering for plotting
row_reordering <-
    get_row_reordering(
        row_clustering = hc_row,
        chain = chains,
        burnin = burnin,
        dat = purrr::map(my_dat, ~ dplyr::select(.x, 4:last_col())),
        multichain = TRUE
    )

molten <- bind_rows(my_dat) %>%
    dplyr::mutate(row = row_reordering) %>%
    dplyr::select(4:last_col()) %>%
    tidyr::pivot_longer(!last_col(), names_to = "cell") %>%
    # Relevel factors, for column sorting on the plot
    mutate(cell = forcats::fct_relevel(cell, ~ hc_col$labels[hc_col$order]))


p1 <- ggplot(data = molten,
             aes(x = cell,
                 y = row,
                 fill = value)) +
  geom_raster() +
  theme_minimal() +
  theme(
    axis.text.x = element_blank(),
  ) +
  labs(fill = "Z-score", x = "", y = "") +
  ggtitle("Bi-clustering heatmap") +
  scale_fill_distiller(palette = "Greens", direction = 1) +
  coord_flip()
print(p1)

#-------------------------------------------------------------------------------
# Read in colormap for plotting, to match cell type by function 
#-------------------------------------------------------------------------------
pal <- read_delim(file = "data/color_mapper.txt", col_names = FALSE, delim = " ") %>%
    set_names(c("cell_pop", "hex"))

# mean trends ---------------------------------------------------
# each row of merged_mu should correspond to a different factor
# each column is a cell population
mu <- merged$merged_mu %>%
    as_tibble(.name_repair = "unique") %>%
    set_names(cell_pops_to_keep) %>%
    mutate(group = 1:n(), .before = 1) %>%
    pivot_longer(cols = 2:last_col(), names_to = "cell_pop", values_to = "mu_est")

# covariance trends ---------------------------------------------------
sigmas <- merged$merged_sigma %>%
    apply(MARGIN = 3, function(X) {
      X %>%
        set_colnames(cell_pops_to_keep) %>%
        as_tibble(X, .name_repair = "minimal")
    })

sigma_df <- sigmas %>%
    imap_dfr(
        ~ mutate(
            .x,
            group = .y,
            cell_pop1 = cell_pops_to_keep,
            .before = 1
        ) %>%
            pivot_longer(
                cols = 3:last_col(),
                names_to = "cell_pop2",
                values_to = "covariance"
            )
    )

#-------------------------------------------------------------------------------
# Use row and column clustering for row reordering
#-------------------------------------------------------------------------------
cl <- list(row_clustering = hc_row, col_clustering = hc_col)
pcs <- map(sigmas, ~ prcomp(.x, center = TRUE))

out_plots <- list()
clusts_to_plot <- seq_along(pcs)

for(ccc in clusts_to_plot) {
    pc_df <- as_tibble(pcs[[ccc]]$rotation) %>%
        mutate(cell_pop = cell_pops_to_keep, .before = 1) %>%
        pivot_longer(cols = !cell_pop, names_to = "PC", values_to = "score") %>%
        mutate(PC = factor(PC, levels = paste0("PC", seq_along(cell_pops_to_keep)))) %>%
        group_split(PC) %>%
        map2_dfr((pcs[[ccc]]$sdev ^2) / sum(pcs[[ccc]]$sdev^2), ~ mutate(.x, percent_var = round(.y * 100, digits = 2))) %>%
        left_join(
            filter(mu, group == ccc) %>%
                mutate(group = NULL), by = "cell_pop") %>%
        left_join(pal, by = "cell_pop") %>%
        mutate(
            cell_pop = factor(cell_pop, levels = cl$col_clustering$labels[cl$col_clustering$order])
        ) %>%
        arrange(cell_pop) %>%
        mutate(hex = factor(hex, levels = unique(.$hex)))
    
   
    mu_plot <- ggplot(filter(pc_df, PC %in% "PC1"), aes(x = cell_pop, y  = mu_est)) +
        geom_bar(aes(fill = hex), stat = "identity", color = "black", show.legend = FALSE) +
        scale_fill_manual(values = levels(pc_df$hex)) +
        labs(x = "", y = "Estimated\ncluster mean") +
        theme_minimal()
    if(ccc == 12) {
        mu_plot <- mu_plot + theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
                                   legend.position = "none")
    } else {
        mu_plot <- mu_plot + theme(axis.text.x = element_blank(), legend.position = "none")
    }
        
        
    pc_df2 <- filter(pc_df, PC %in% paste0("PC", 1:2)) %>%
        unite(col = PC_var, PC, percent_var, sep = " (") %>%
        mutate(PC_var = paste0(PC_var, "%)"))
    
    pc_plot <- ggplot(pc_df2, aes(x = cell_pop, y = score)) +
        geom_bar(aes(fill = hex), stat = "identity", color = "black", show.legend = FALSE) +
        scale_fill_manual(values = levels(pc_df$hex)) +
        labs(x = "") +
        facet_wrap(~ PC_var, nrow = 1, ncol = 2) +
        theme_minimal()
    if(ccc == 12) {
        pc_plot <- pc_plot + theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
                                   legend.position = "none")
    } else {
        pc_plot <- pc_plot + theme(axis.text.x = element_blank(), legend.position = "none")
    }
    
    out_plots[[ccc]] <- cowplot::plot_grid(mu_plot, pc_plot, nrow = 1, ncol = 2, rel_widths = c(1,2))
}

cowplot::plot_grid(plotlist = out_plots, nrow = length(out_plots), ncol = 1, rel_heights = c(rep(1,11), 2))