R/homer_heatmap.R

In TedCCLeung/PhotoperiodMotif: What the Package Does (One Line, Title Case)

input_homer <- function(
  dir = "/Users/TedCCLeung/Documents/Projects/Photoperiod/2_analysis/2_pipeline/PhotoperiodMotif/homer_08/"
){
  motif_files <- list.files(path = dir, recursive = TRUE, pattern = "knownResults.txt", all.files = TRUE, full.names = TRUE)
  cond_names <- motif_files %>% strsplit(split = "/") %>% sapply(function(x){magrittr::extract2(x, length(x)-1)})

  df_result <- Reduce(rbind, Map(read_homer_output, motif_files, cond_names))

  return(df_result)
}

process_homer_matrix <- function(
  dir,
  df_result,
  motifs = universalmotif::read_homer("/Users/TedCCLeung/Documents/Projects/Packages/all.8.motif"),
  pval_threshold = 1e-3,
  cutree_height = 0.5
){

  ## Make the matrices
  mat_pval <- df_result[, c("name", "sample", "pvalue")] %>%
    tidyr::pivot_wider(names_from = "sample", values_from = "pvalue") %>%
    tibble::column_to_rownames(var = "name") %>%
    log10() %>%
    as.matrix() %>%
    dplyr::na_if(Inf) %>%
    dplyr::na_if(-Inf) %>%
    tidyr::replace_na(0)

  mat_fdr <- df_result[, c("name", "sample", "fdr")] %>%
    tidyr::pivot_wider(names_from = "sample", values_from = "fdr") %>%
    tibble::column_to_rownames(var = "name") %>%
    log10() %>%
    as.matrix() %>%
    dplyr::na_if(Inf) %>%
    dplyr::na_if(-Inf) %>%
    tidyr::replace_na(0)

  mat_fc <- df_result[, c("name", "sample", "fold_change")] %>%
    tidyr::pivot_wider(names_from = "sample", values_from = "fold_change") %>%
    tibble::column_to_rownames(var = "name") %>%
    log2() %>%
    suppressWarnings() %>%
    as.matrix() %>%
    dplyr::na_if(Inf) %>%
    dplyr::na_if(-Inf) %>%
    tidyr::replace_na(0)

  ## Filter
  row_min <- apply(mat_pval, 1, min)
  common_motifs <- row_min[row_min < log10(pval_threshold)] %>% names
  motifs <- universalmotif::filter_motifs(motifs, name = common_motifs)

  mat_pval <- mat_pval[, colnames(mat_pval) %in% common_motifs]
  mat_fdr <- mat_fdr[, colnames(mat_fdr) %in% common_motifs]
  mat_fc <- mat_fc[, colnames(mat_fc) %in% common_motifs]

  ## Clustering ---
  hclust_out <- motifs %>% hclust_motifs() %>% suppressWarnings()
  cutree_out <- stats::cutree(hclust_out, h = cutree_height)

  ## Heat maps ---
  heatmap_PVAL <- heatmap_from_Fisher(mat_pval, hclust_out, cutree_out, motifs, "PVAL")
  heatmap_FDR <- heatmap_from_Fisher(mat_fdr, hclust_out, cutree_out, motifs, "FDR")
  heatmap_FC <- heatmap_from_Fisher(mat_fc, hclust_out, cutree_out, motifs, "FC")

  ## Get the motif families for plotting the bar ----
  ## Any of the heatmap can be used for the main
  heatmap_TF_family <- heatmap_motif_to_CISBP_family(mat_fc, hclust_out, motifs, main = heatmap_FC) %>% heatmap_TF_family()
  heatmap_motif_source <- heatmap_motif_source(matrix = mat_fc, hclust_out = hclust_out, main = heatmap_FC)

  ## Add the extra information to the heat maps ----
  heatmap_PVAL <- heatmap_PVAL + heatmap_TF_family + heatmap_motif_source
  heatmap_FDR <- heatmap_FDR + heatmap_TF_family + heatmap_motif_source
  heatmap_FC <- heatmap_FC + heatmap_TF_family + heatmap_motif_source

  ## Output pdf ----
  return(list(heatmap_PVAL = heatmap_PVAL, heatmap_FDR = heatmap_FDR, heatmap_FC = heatmap_FC))
}