plotHeatmap_trial

workflowr

• Summary
• Checks
• Past versions

Last updated: 2026-03-20

Checks: 7 0

Knit directory: DXR_continue/

This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

---

R Markdown file: up-to-date

Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20250701)

The command set.seed(20250701) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: f8dc9bf

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version f8dc9bf. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:

Ignored files:
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    analysis/figure/
    Ignored:    data/Bed_exports/
    Ignored:    data/Cormotif_data/
    Ignored:    data/DER_data/
    Ignored:    data/Other_paper_data/
    Ignored:    data/RDS_files/
    Ignored:    data/TE_annotation/
    Ignored:    data/alignment_summary.txt
    Ignored:    data/all_peak_final_dataframe.txt
    Ignored:    data/cell_line_info_.tsv
    Ignored:    data/full_summary_QC_metrics.txt
    Ignored:    data/motif_lists/
    Ignored:    data/number_frag_peaks_summary.txt

Untracked files:
    Untracked:  H3K27ac_all_regions_test.bed
    Untracked:  H3K27ac_consensus_clusters_test.bed
    Untracked:  analysis/GREAT_H3K27ac.Rmd
    Untracked:  analysis/H3K27ac_ChromHMM_FC.Rmd
    Untracked:  analysis/H3K27ac_TEtype.Rmd
    Untracked:  analysis/H3K27me3_TE_investigation.Rmd
    Untracked:  analysis/H3K36me3_TE_investigation.Rmd
    Untracked:  analysis/Top2a_Top2b_expression.Rmd
    Untracked:  analysis/maps_and_plots.Rmd
    Untracked:  analysis/proteomics.Rmd
    Untracked:  code/For_john.R
    Untracked:  other_analysis/

Unstaged changes:
    Modified:   analysis/H3K27_TE_overlap.Rmd
    Modified:   analysis/H3K27ac_cisRE.Rmd
    Modified:   analysis/H3K27ac_summit_processing.Rmd
    Modified:   analysis/H3K27me3_summit_processing.Rmd
    Modified:   analysis/H3K36me3_summit_processing.Rmd
    Modified:   analysis/H3K9me3_TE_investigation.Rmd
    Modified:   analysis/chromHMM.Rmd
    Modified:   analysis/dual_histone_TE_investigation.Rmd
    Modified:   analysis/summit_files_processing.Rmd

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

---

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/computeMatrixplots.Rmd) and HTML (docs/computeMatrixplots.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File	Version	Author	Date	Message
Rmd	f8dc9bf	reneeisnowhere	2026-03-20	updates
html	8fe7ed0	reneeisnowhere	2026-03-19	Build site.
Rmd	e77cb37	reneeisnowhere	2026-03-19	first commit

---

library(tidyverse)
library(GenomicRanges)
library(plyranges)
library(genomation)
library(readr)
library(rtracklayer)
library(stringr)
library(ggrepel)
library(DT)
library(data.table)
library(circlize)
library(ComplexHeatmap)

repeatmasker <- read_delim("data/Other_paper_data/repeatmasker_20250911.txt", 
    delim = "\t", escape_double = FALSE, 
    trim_ws = TRUE)
colnames(repeatmasker)

 [1] "#bin"      "swScore"   "milliDiv"  "milliDel"  "milliIns"  "genoName" 
 [7] "genoStart" "genoEnd"   "genoLeft"  "strand"    "repName"   "repClass" 
[13] "repFamily" "repStart"  "repEnd"    "repLeft"   "id"       

autosomes <- paste0("chr", 1:22)
repeatmasker_clean <- repeatmasker %>% mutate(
    strand = ifelse(strand == "C", "-", "+")
  ) %>% 
   mutate(
    start = genoStart + 1,
    end   = genoEnd)%>% 
  mutate(repFamily= str_remove(repFamily, "\\?$")) %>% 
  dplyr::filter(genoName %in% autosomes) %>% 
  mutate(RM_id=paste0(genoName,":",start,"-",end,":",id))


rpt_split <- split(repeatmasker_clean, repeatmasker_clean$repClass)

rpt_split_gr_list <- lapply(rpt_split, function(df) {
  GRanges(
    seqnames = df$genoName,
    ranges   = IRanges(start = df$start, end = df$end),
    strand   = df$strand,
    repName   = df$repName,
    repClass  = df$repClass,
    repFamily = df$repFamily,
    swScore   = df$swScore,
    milliDiv  = df$milliDiv,
    milliDel  = df$milliDel,
    milliIns  = df$milliIns,
    RM_id        = df$RM_id
  )
})

matrix_file <- "C:/Users/renee/Other_projects_data/DXR_data/computeMatrix_output/four_histone/matrix_MER61E_ERV1_LTR.gz"
cm <- fread(cmd = paste("gzip -dc", matrix_file), skip = 1)
dim(cm)

[1]   261 21606

head(cm[,1:10])

       V1       V2       V3                       V4    V5     V6       V7
   <char>    <int>    <int>                   <char> <num> <char>    <num>
1:   chr1 63032707 63032840 chr1:63032708-63032840:1     0      + 0.000000
2:   chr1 63033079 63033512 chr1:63033080-63033512:1     0      + 0.000000
3:   chr1 69055251 69055817 chr1:69055252-69055817:1     0      + 0.000000
4:   chr1 71345730 71345939 chr1:71345731-71345939:1     0      + 0.000000
5:   chr1 72000824 72001353 chr1:72000825-72001353:1     0      + 0.000000
6:   chr1 76040192 76040658 chr1:76040193-76040658:1     0      + 0.014045
         V8       V9      V10
      <num>    <num>    <num>
1: 0.000000 0.000000 0.000000
2: 0.000000 0.025281 0.028090
3: 0.000000 0.000000 0.000000
4: 0.000000 0.000000 0.000000
5: 0.000000 0.000000 0.000000
6: 0.014045 0.014045 0.014045

sample_names <- c(
"H3K27ac_DOX_144R_merged","H3K27ac_DOX_24R_merged","H3K27ac_DOX_24T_merged",
"H3K27ac_VEH_144R_merged","H3K27ac_VEH_24R_merged","H3K27ac_VEH_24T_merged",
"H3K27me3_DOX_144R_merged","H3K27me3_DOX_24R_merged","H3K27me3_DOX_24T_merged",
"H3K27me3_VEH_144R_merged","H3K27me3_VEH_24R_merged","H3K27me3_VEH_24T_merged",
"H3K36me3_DOX_144R_merged","H3K36me3_DOX_24R_merged","H3K36me3_DOX_24T_merged",
"H3K36me3_VEH_144R_merged","H3K36me3_VEH_24R_merged","H3K36me3_VEH_24T_merged",
"H3K9me3_DOX_144R_merged","H3K9me3_DOX_24R_merged","H3K9me3_DOX_24T_merged",
"H3K9me3_VEH_144R_merged","H3K9me3_VEH_24R_merged","H3K9me3_VEH_24T_merged"
)
meta <- meta <- cm[,1:6]
bins <- 900
signal_list <- lapply(1:length(sample_names), function(i) {
  start <- 6 + (i-1)*bins + 1
  end <- 6 + i*bins
  as.matrix(cm[,start:end])
})

names(signal_list) <- sample_names

profiles <- lapply(signal_list, colMeans)

profiles_df <- do.call(cbind, profiles)
profiles_df <- as.data.frame(profiles_df)

profiles_df$bin <- 1:nrow(profiles_df)

long <- pivot_longer(
  profiles_df,
  -bin,
  names_to="sample",
  values_to="signal"
) %>%
  separate(sample, into=c("mark","condition","time"), sep="_", remove = FALSE)


ggplot(long, aes(bin, signal, color=sample)) +
  geom_line(size=1) +
  theme_bw() +
  labs(
    x="Scaled TE region",
    y="ChIP signal"
  )+
facet_wrap(~time)

Past versions of unnamed-chunk-3-1.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

long$condition <- ifelse(grepl("DOX", long$sample),"DOX","VEH")
long$mark <- sub("_.*","",long$sample)


ggplot(long, aes(bin, signal, color=time)) +
  geom_line(size=1) +
  facet_grid(mark ~ condition) +
  theme_bw()+
  geom_vline(xintercept=c(200,700), linetype="dashed")+
  ggtitle("signal across all LTR:ERV1:MER61E")

Past versions of unnamed-chunk-3-2.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

plot_deeptools_matrix <- function(matrix_file, bins = 900){

  library(data.table)
  library(dplyr)
  library(tidyr)
  library(ggplot2)
  library(stringr)

  message("Processing: ", matrix_file)

  cm <- fread(cmd = paste("gzip -dc", matrix_file), skip = 1)

  sample_names <- c(
  "H3K27ac_DOX_144R_merged","H3K27ac_DOX_24R_merged","H3K27ac_DOX_24T_merged",
  "H3K27ac_VEH_144R_merged","H3K27ac_VEH_24R_merged","H3K27ac_VEH_24T_merged",
  "H3K27me3_DOX_144R_merged","H3K27me3_DOX_24R_merged","H3K27me3_DOX_24T_merged",
  "H3K27me3_VEH_144R_merged","H3K27me3_VEH_24R_merged","H3K27me3_VEH_24T_merged",
  "H3K36me3_DOX_144R_merged","H3K36me3_DOX_24R_merged","H3K36me3_DOX_24T_merged",
  "H3K36me3_VEH_144R_merged","H3K36me3_VEH_24R_merged","H3K36me3_VEH_24T_merged",
  "H3K9me3_DOX_144R_merged","H3K9me3_DOX_24R_merged","H3K9me3_DOX_24T_merged",
  "H3K9me3_VEH_144R_merged","H3K9me3_VEH_24R_merged","H3K9me3_VEH_24T_merged"
  )

  signal_list <- lapply(seq_along(sample_names), function(i){

    start <- 6 + (i-1)*bins + 1
    end <- 6 + i*bins

    as.matrix(cm[,start:end])
  })

  names(signal_list) <- sample_names

  profiles <- lapply(signal_list, colMeans)

  profiles_df <- do.call(cbind, profiles) |> as.data.frame()
  profiles_df$bin <- 1:nrow(profiles_df)

  long <- profiles_df |>
    pivot_longer(-bin, names_to="sample", values_to="signal") |>
    separate(sample, into=c("mark","condition","time","extra"), sep="_") |>
    mutate(sample = paste(mark,condition,time,sep="_"))

  # remove unwanted marks
  long <- long |> filter(!mark %in% c("H3K27me3","H3K36me3"))

  title <- tools::file_path_sans_ext(basename(matrix_file))

  p <- ggplot(long, aes(bin, signal, color=time)) +
    geom_line(linewidth=1) +
    facet_grid(mark ~ condition) +
    geom_vline(xintercept=c(200,700), linetype="dashed") +
    theme_bw() +
    ggtitle(title)

  return(p)
}

plot_deeptools_matrix("C:/Users/renee/Other_projects_data/DXR_data/computeMatrix_output/four_histone/matrix_MER61E_ERV1_LTR.gz")

Past versions of unnamed-chunk-5-1.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

matrix_dir <- "C:/Users/renee/Other_projects_data/DXR_data/computeMatrix_output/four_histone/"

matrix_files <- list.files(
  matrix_dir,
  pattern="matrix_.*gz$",
  full.names=TRUE
)

plots <- lapply(matrix_files, plot_deeptools_matrix)

plots

[[1]]

Past versions of unnamed-chunk-7-1.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

[[2]]

Past versions of unnamed-chunk-7-2.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

[[3]]

Past versions of unnamed-chunk-7-3.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

[[4]]

Past versions of unnamed-chunk-7-4.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

[[5]]

Past versions of unnamed-chunk-7-5.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

[[6]]

Past versions of unnamed-chunk-7-6.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

[[7]]

Past versions of unnamed-chunk-7-7.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

[[8]]

Past versions of unnamed-chunk-7-8.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

[[9]]

Past versions of unnamed-chunk-7-9.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

[[10]]

Past versions of unnamed-chunk-7-10.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

[[11]]

Past versions of unnamed-chunk-7-11.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

[[12]]

Past versions of unnamed-chunk-7-12.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

[[13]]

Past versions of unnamed-chunk-7-13.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

[[14]]

Past versions of unnamed-chunk-7-14.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

[[15]]

Past versions of unnamed-chunk-7-15.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

[[16]]

Past versions of unnamed-chunk-7-16.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

[[17]]

Past versions of unnamed-chunk-7-17.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

[[18]]

Past versions of unnamed-chunk-7-18.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

[[19]]

Past versions of unnamed-chunk-7-19.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

[[20]]

Past versions of unnamed-chunk-7-20.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Heatmap plotting

meta <- cm[,1:6]
signal <- cm[,7:ncol(cm)]

bins <- 900
sample_index <- 2

start <- (sample_index - 1) * bins + 1
end <- sample_index * bins

mat <- as.matrix(signal[,start:end])
### sorting rows by signal
row_order <- order(rowMeans(mat), decreasing = TRUE)
mat <- mat[row_order,]


### Color scale
col_fun <- colorRamp2(
  c(0, 0.5, 2),
  c("white", "orange", "red")
)
#### color scale for centered data
# col_fun <- colorRamp2(
#   c(-2, 0, 2),
#   c("blue", "white", "red")
# )

region <- c(
  rep("Upstream",200),
  rep("TE",500),
  rep("Downstream",200)
)
Heatmap(
  mat,
  col = col_fun,
  name = "signal",
  show_row_names = FALSE,
  show_column_names = FALSE,
  column_split = region,
  cluster_columns = FALSE
)

Past versions of unnamed-chunk-8-1.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

plot_2hist_deeptools_matrix <- function(matrix_file, bins = 400){

  library(data.table)
  library(dplyr)
  library(tidyr)
  library(ggplot2)
  library(stringr)

  message("Processing: ", matrix_file)

  cm <- fread(cmd = paste("gzip -dc", matrix_file), skip = 1)

  sample_names <- c(
  "H3K27ac_DOX_144R_merged","H3K27ac_DOX_24R_merged","H3K27ac_DOX_24T_merged",
  "H3K27ac_VEH_144R_merged","H3K27ac_VEH_24R_merged","H3K27ac_VEH_24T_merged",
  "H3K9me3_DOX_144R_merged","H3K9me3_DOX_24R_merged","H3K9me3_DOX_24T_merged",
  "H3K9me3_VEH_144R_merged","H3K9me3_VEH_24R_merged","H3K9me3_VEH_24T_merged"
  )

  signal_list <- lapply(seq_along(sample_names), function(i){

    start <- 6 + (i-1)*bins + 1
    end <- 6 + i*bins

    as.matrix(cm[,start:end])
  })

  names(signal_list) <- sample_names

  profiles <- lapply(signal_list, colMeans)

  profiles_df <- do.call(cbind, profiles) |> as.data.frame()
 profiles_df$pos_kb <- seq(-2, 2, length.out = nrow(profiles_df))
 
 
  long <- profiles_df |>
    pivot_longer(-pos_kb, names_to="sample", values_to="signal") |>
    separate(sample, into=c("mark","condition","time","extra"), sep="_") |>
    mutate(sample = paste(mark,condition,time,sep="_"))

  # # remove unwanted marks
  # long <- long |> filter(!mark %in% c("H3K27me3","H3K36me3"))

  title <- tools::file_path_sans_ext(basename(matrix_file))

  p <- ggplot(long, aes(pos_kb, signal, color=time)) +
    geom_line(linewidth=1) +
    facet_grid(mark ~ condition) +
    geom_vline(xintercept = 0, linetype="dashed")+
    # geom_smooth(se=FALSE, method = "loess")+
    theme_bw() +
    ggtitle(title)

  return(p)
}

plot_2hist_deeptools_matrix("C:/Users/renee/Other_projects_data/DXR_data/computeMatrix_output_centered2histone/matcenter_MER61E_ERV1_LTR.gz", bins=400)

Past versions of unnamed-chunk-9-1.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

matrix_dir_cent <- "C:/Users/renee/Other_projects_data/DXR_data/computeMatrix_output_centered2histone/"

matrix_files_cent <- list.files(
  matrix_dir_cent,
  pattern="matcenter_.*gz$",
  full.names=TRUE
)

plots_cent <- lapply(matrix_files_cent, plot_2hist_deeptools_matrix)
for (p in plots_cent) {
  print(p)
}

Past versions of unnamed-chunk-10-1.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-10-2.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-10-3.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-10-4.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-10-5.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-10-6.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-10-7.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-10-8.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-10-9.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-10-10.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-10-11.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-10-12.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-10-13.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-10-14.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-10-15.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-10-16.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-10-17.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-10-18.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-10-19.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

col_fun <- colorRamp2(
  c(-2, 0, 2),
  c("blue", "white", "red")
)
bins=400
cm <- fread(cmd = paste("gzip -dc", "C:/Users/renee/Other_projects_data/DXR_data/computeMatrix_output_centered2histone/matcenter_MER61E_ERV1_LTR.gz"), skip = 1)

  sample_names <- c(
  "H3K27ac_DOX_144R_merged","H3K27ac_DOX_24R_merged","H3K27ac_DOX_24T_merged",
  "H3K27ac_VEH_144R_merged","H3K27ac_VEH_24R_merged","H3K27ac_VEH_24T_merged",
  "H3K9me3_DOX_144R_merged","H3K9me3_DOX_24R_merged","H3K9me3_DOX_24T_merged",
  "H3K9me3_VEH_144R_merged","H3K9me3_VEH_24R_merged","H3K9me3_VEH_24T_merged"
  )
  ##extract per-sample matrices
  signal_list <- lapply(seq_along(sample_names), function(i){

    start <- 6 + (i-1)*bins + 1
    end <- 6 + i*bins

    as.matrix(cm[,start:end])
  })
ref_mat <-signal_list[[3]] 
row_order <- order(rowMeans(ref_mat), decreasing = TRUE)
  names(signal_list) <- sample_names

  for (i in seq_along(signal_list)) {  
  #scaling
  mat <- signal_list[[i]]
mat_scaled <- t(scale(t(mat)))
### sorting rows by signal
# row_order <- order(rowMeans(mat), decreasing = TRUE)
mat_scaled <- mat_scaled[row_order,]

mat_scaled[is.na(mat_scaled)] <- 0

region <- ifelse(seq_len(ncol(mat_scaled)) <= ncol(mat_scaled)/2, "Upstream", "Downstream")


Heatmap(
  mat_scaled,
      col = col_fun,
      name = "signal",
      show_row_names = FALSE,
      show_column_names = FALSE,
      cluster_rows = FALSE,           # Important to preserve row_order
      cluster_columns = FALSE,
      column_split = region,
      column_title = names(signal_list)[i])
}

plot_deeptools_heatmaps <- function(matrix_file,
                                    sample_order_H3K27ac = c("H3K27ac_DOX_24T_merged",
                                                              "H3K27ac_DOX_24R_merged",
                                                              "H3K27ac_DOX_144R_merged",
                                                              "H3K27ac_VEH_24T_merged",
                                                              "H3K27ac_VEH_24R_merged",
                                                              "H3K27ac_VEH_144R_merged"),
                                    sample_order_H3K9me3 = c("H3K9me3_DOX_24T_merged",
                                                              "H3K9me3_DOX_24R_merged",
                                                              "H3K9me3_DOX_144R_merged",
                                                              "H3K9me3_VEH_24T_merged",
                                                              "H3K9me3_VEH_24R_merged",
                                                              "H3K9me3_VEH_144R_merged"),
                                    bins = 400,
                                    reference_sample = "H3K27ac_DOX_24T_merged") {
  
  # Read the matrix
  cm <- fread(cmd = paste("gzip -dc", matrix_file), skip = 1)
  
  # Full list of samples
  sample_names <- c(
    "H3K27ac_DOX_144R_merged","H3K27ac_DOX_24R_merged","H3K27ac_DOX_24T_merged",
    "H3K27ac_VEH_144R_merged","H3K27ac_VEH_24R_merged","H3K27ac_VEH_24T_merged",
    "H3K9me3_DOX_144R_merged","H3K9me3_DOX_24R_merged","H3K9me3_DOX_24T_merged",
    "H3K9me3_VEH_144R_merged","H3K9me3_VEH_24R_merged","H3K9me3_VEH_24T_merged"
  )
  
  # Extract per-sample matrices
  signal_list <- lapply(seq_along(sample_names), function(i){
    start <- 6 + (i-1)*bins + 1
    end <- 6 + i*bins
    as.matrix(cm[, start:end])
  })
  names(signal_list) <- sample_names
  file_part <- str_extract(basename(matrix_file), "(?<=matcenter_).*?(?=_ERV1_LTR)")
  # Determine row order from reference sample
  ref_mat <- signal_list[[which(names(signal_list) == reference_sample)]]
  row_order <- order(rowMeans(ref_mat), decreasing = TRUE)
  
  # Color function
  col_fun <- colorRamp2(c(-2, 0, 2), c("blue", "white", "red"))
  
  ## Helper function to scale and reorder matrix
  prepare_mat <- function(mat) {
    mat_scaled <- t(scale(t(mat)))          # row scaling
    mat_scaled[is.na(mat_scaled)] <- 0
    mat_scaled[row_order, , drop = FALSE]   # apply reference row order
  }
  
  # Function to create heatmaps for a given sample vector
 plot_heatmap_group <- function(sample_vector, group_name) {
    ht_list <- NULL
    
    for (s in sample_vector) {
      mat <- prepare_mat(signal_list[[s]])
      region <- ifelse(seq_len(ncol(mat)) <= ncol(mat)/2, "Upstream", "Downstream")
      
      # Shorten the column title
      short_title <- str_extract(s, "(DOX|VEH)_\\d+[RT]")        # keeps "DOX_24T" etc.
      short_title <- str_replace(short_title, "_merged", "")
      
      
      ht <- Heatmap(mat,
                    col = col_fun,
                    name = "signal",
                    show_row_names = FALSE,
                    show_column_names = FALSE,
                    cluster_rows = FALSE,
                    cluster_columns = FALSE,
                    # column_split = region,
                    column_title = short_title)
      
      if (is.null(ht_list)) {
        ht_list <- ht
      } else {
        ht_list <- ht_list + ht
      }
    }
    
    # Draw all heatmaps with a master title
    draw(ht_list, heatmap_legend_side = "right", column_title = paste0(group_name," - ",file_part))
  }
  
  # Plot H3K27ac group
  plot_heatmap_group(sample_order_H3K27ac, "H3K27ac")
  
  # Plot H3K9me3 group
  plot_heatmap_group(sample_order_H3K9me3, "H3K9me3")
}

plot_deeptools_heatmaps("C:/Users/renee/Other_projects_data/DXR_data/computeMatrix_output_centered2histone/matcenter_MER61E_ERV1_LTR.gz")

Past versions of unnamed-chunk-13-1.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-2.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

for (f in matrix_files_cent) {
  plot_deeptools_heatmaps(f)   # plots automatically
}

Past versions of unnamed-chunk-13-3.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-4.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-5.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-6.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-7.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-8.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-9.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-10.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-11.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-12.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-13.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-14.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-15.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-16.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-17.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-18.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-19.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-20.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-21.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-22.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-23.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-24.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-25.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-26.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-27.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-28.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-29.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-30.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-31.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-32.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-33.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-34.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-35.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-36.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-37.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-38.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-39.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

Past versions of unnamed-chunk-13-40.png

Version	Author	Date
8fe7ed0	reneeisnowhere	2026-03-19

finding switches

plot_switching_TEs <- function(matrix_file,
                               bins = 400,
                               fc_threshold_ac = 0.01,
                               fc_threshold_me3 = -0.5) {
  
  # -----------------------------
  # 1️⃣ Load matrix
  # -----------------------------
  cm <- fread(cmd = paste("gzip -dc", matrix_file), skip = 1)
  
  sample_names <- c(
    "H3K27ac_DOX_144R_merged","H3K27ac_DOX_24R_merged","H3K27ac_DOX_24T_merged",
    "H3K27ac_VEH_144R_merged","H3K27ac_VEH_24R_merged","H3K27ac_VEH_24T_merged",
    "H3K9me3_DOX_144R_merged","H3K9me3_DOX_24R_merged","H3K9me3_DOX_24T_merged",
    "H3K9me3_VEH_144R_merged","H3K9me3_VEH_24R_merged","H3K9me3_VEH_24T_merged"
  )
  
  # -----------------------------
  # 2️⃣ Extract matrices
  # -----------------------------
  signal_list <- lapply(seq_along(sample_names), function(i){
    start <- 6 + (i-1)*bins + 1
    end <- 6 + i*bins
    as.matrix(cm[, start:end])
  })
  names(signal_list) <- sample_names
  
  # -----------------------------
  # 3️⃣ Compute average signal across bins for each mark & condition
  # -----------------------------
  avg_signal <- function(mark, cond) {
    rows <- grep(paste0(mark,"_",cond), names(signal_list))
    # Average across replicates
    Reduce("+", signal_list[rows]) / length(rows)
  }
  
  H3K27ac_DOX_avg <- avg_signal("H3K27ac", "DOX")
  H3K27ac_VEH_avg <- avg_signal("H3K27ac", "VEH")
  H3K9me3_DOX_avg <- avg_signal("H3K9me3", "DOX")
  H3K9me3_VEH_avg <- avg_signal("H3K9me3", "VEH")
  
  # -----------------------------
  # 4️⃣ Compute log2 fold-changes
  # -----------------------------
  H3K27ac_log2fc <- log2((H3K27ac_DOX_avg + 1) / (H3K27ac_VEH_avg + 1))
  H3K9me3_log2fc <- log2((H3K9me3_DOX_avg + 1) / (H3K9me3_VEH_avg + 1))
  
  # -----------------------------
  # 5️⃣ Identify switching TEs
  # -----------------------------
  switching_rows <- which(H3K27ac_log2fc > fc_threshold_ac &
                          H3K9me3_log2fc < fc_threshold_me3)
  
  if(length(switching_rows) == 0) {
    message("No switching TEs found with these thresholds.")
    return(NULL)
  }
  
  # -----------------------------
  # 6️⃣ Prepare matrices for plotting
  # -----------------------------
  prepare_mat <- function(mat, rows) {
    mat_sub <- mat[rows, , drop = FALSE]
    mat_scaled <- t(scale(t(mat_sub)))
    mat_scaled[is.na(mat_scaled)] <- 0
    mat_scaled
  }
  
  # Reference row order: H3K27ac DOX 24T
  ref_mat <- signal_list[["H3K27ac_DOX_24T_merged"]][switching_rows, ]
  row_order <- order(rowMeans(ref_mat), decreasing = TRUE)
  
  # Function to make heatmaps for a given mark
  make_heatmap <- function(mark, samples) {
    ht_list <- NULL
    for(s in samples) {
      mat_scaled <- prepare_mat(signal_list[[s]], switching_rows)
      mat_scaled <- mat_scaled[row_order, ]
      short_title <- str_extract(s, "(DOX|VEH)_\\d+[RT]")
      ht <- Heatmap(mat_scaled,
                    col = colorRamp2(c(-2,0,2), c("blue","white","red")),
                    name = "signal",
                    show_row_names = FALSE,
                    show_column_names = FALSE,
                    cluster_rows = FALSE,
                    cluster_columns = FALSE,
                    column_title = short_title)
      ht_list <- if(is.null(ht_list)) ht else ht_list + ht
    }
    ht_list
  }
  
  # Define sample order for plotting
  H3K27ac_samples <- c("H3K27ac_DOX_24T_merged","H3K27ac_DOX_24R_merged","H3K27ac_DOX_144R_merged",
                       "H3K27ac_VEH_24T_merged","H3K27ac_VEH_24R_merged","H3K27ac_VEH_144R_merged")
  H3K9me3_samples <- c("H3K9me3_DOX_24T_merged","H3K9me3_DOX_24R_merged","H3K9me3_DOX_144R_merged",
                       "H3K9me3_VEH_24T_merged","H3K9me3_VEH_24R_merged","H3K9me3_VEH_144R_merged")
  
  ht_H3K27ac <- make_heatmap("H3K27ac", H3K27ac_samples)
  ht_H3K9me3 <- make_heatmap("H3K9me3", H3K9me3_samples)
  
  # -----------------------------
  # 7️⃣ Extract file part for master title
  # -----------------------------
  file_part <- str_extract(basename(matrix_file), "(?<=matcenter_).*?(?=_ERV1_LTR)")
  
  draw(ht_H3K27ac, heatmap_legend_side = "right", column_title = paste0("H3K27ac - ", file_part))
  draw(ht_H3K9me3, heatmap_legend_side = "right", column_title = paste0("H3K9me3 - ", file_part))
  
  return(list(H3K27ac = ht_H3K27ac, H3K9me3 = ht_H3K9me3, switching_rows = switching_rows))
}

result <- plot_switching_TEs(
  "C:/Users/renee/Other_projects_data/DXR_data/computeMatrix_output_centered2histone/matcenter_MER61E_ERV1_LTR.gz")

# 
##debug
matrix_file <- ("C:/Users/renee/Other_projects_data/DXR_data/computeMatrix_output_centered2histone/matcenter_MER61E_ERV1_LTR.gz")

bins <- 400

H3K27ac_DOX_24T <- rowMeans(signal_list[["H3K27ac_DOX_24T_merged"]])
H3K27ac_VEH_24T <- rowMeans(signal_list[["H3K27ac_VEH_24T_merged"]])

H3K9me3_DOX_24T <- rowMeans(signal_list[["H3K9me3_DOX_24T_merged"]])
H3K9me3_VEH_24T <- rowMeans(signal_list[["H3K9me3_VEH_24T_merged"]])


summary(H3K27ac_DOX_24T)
summary(H3K27ac_VEH_24T)

summary(H3K9me3_DOX_24T)
summary(H3K9me3_H3K9me3_VEH_24TVEH_144R)

H3K27ac_log2fc_24T <- log2((H3K27ac_DOX_24T + 1) / (H3K27ac_VEH_24T + 1))
H3K9me3_log2fc_24T <- log2((H3K9me3_DOX_24T + 1) / (H3K9me3_VEH_24T + 1))

summary(H3K27ac_log2fc_24T)
summary(H3K9me3_log2fc_24T)

fc_threshold_ac <- 0.01  
fc_threshold_me3 <- 0    

switching_rows <- which(
  H3K27ac_log2fc_24T > fc_threshold_ac &
  H3K9me3_log2fc_24T <= fc_threshold_me3)


mat_H3K27ac_24T <- signal_list[["H3K27ac_DOX_24T_merged"]][switching_rows, ]
mat_H3K9me3_24T <- signal_list[["H3K9me3_DOX_24T_merged"]][switching_rows, ]

row_order <- order(rowMeans(mat_H3K27ac_24T), decreasing = TRUE)
mat_H3K27ac_24T <- mat_H3K27ac_24T[row_order, ]
mat_H3K9me3_24T <- mat_H3K9me3_24T[row_order, ]

bed_names <- cm[switching_rows, 4]  # column 4 = name
test <- rpt_split$LTR

test %>% 
  dplyr::filter(RM_id %in% bed_names$V4)
H3K27ac_summit_ols %>% dplyr::filter(ID %in%bed_names$V4)

Session information

sessionInfo()

R version 4.4.2 (2024-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26200)

Matrix products: default


locale:
[1] LC_COLLATE=English_United States.utf8 
[2] LC_CTYPE=English_United States.utf8   
[3] LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.utf8    

time zone: America/Chicago
tzcode source: internal

attached base packages:
[1] grid      stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] ComplexHeatmap_2.22.0 circlize_0.4.16       data.table_1.17.8    
 [4] DT_0.33               ggrepel_0.9.6         rtracklayer_1.66.0   
 [7] genomation_1.38.0     plyranges_1.26.0      GenomicRanges_1.58.0 
[10] GenomeInfoDb_1.42.3   IRanges_2.40.1        S4Vectors_0.44.0     
[13] BiocGenerics_0.52.0   lubridate_1.9.4       forcats_1.0.0        
[16] stringr_1.5.1         dplyr_1.1.4           purrr_1.1.0          
[19] readr_2.1.5           tidyr_1.3.1           tibble_3.3.0         
[22] ggplot2_3.5.2         tidyverse_2.0.0       workflowr_1.7.1      

loaded via a namespace (and not attached):
  [1] bitops_1.0-9                rlang_1.1.6                
  [3] magrittr_2.0.3              clue_0.3-66                
  [5] GetoptLong_1.0.5            git2r_0.36.2               
  [7] gridBase_0.4-7              matrixStats_1.5.0          
  [9] compiler_4.4.2              getPass_0.2-4              
 [11] png_0.1-8                   callr_3.7.6                
 [13] vctrs_0.6.5                 reshape2_1.4.4             
 [15] shape_1.4.6.1               pkgconfig_2.0.3            
 [17] crayon_1.5.3                fastmap_1.2.0              
 [19] magick_2.8.7                XVector_0.46.0             
 [21] labeling_0.4.3              Rsamtools_2.22.0           
 [23] promises_1.3.3              rmarkdown_2.29             
 [25] tzdb_0.5.0                  UCSC.utils_1.2.0           
 [27] ps_1.9.1                    bit_4.6.0                  
 [29] xfun_0.52                   zlibbioc_1.52.0            
 [31] cachem_1.1.0                jsonlite_2.0.0             
 [33] later_1.4.2                 DelayedArray_0.32.0        
 [35] BiocParallel_1.40.2         cluster_2.1.8.1            
 [37] parallel_4.4.2              R6_2.6.1                   
 [39] bslib_0.9.0                 stringi_1.8.7              
 [41] RColorBrewer_1.1-3          jquerylib_0.1.4            
 [43] iterators_1.0.14            Rcpp_1.1.0                 
 [45] SummarizedExperiment_1.36.0 knitr_1.50                 
 [47] httpuv_1.6.16               Matrix_1.7-3               
 [49] timechange_0.3.0            tidyselect_1.2.1           
 [51] rstudioapi_0.17.1           dichromat_2.0-0.1          
 [53] abind_1.4-8                 yaml_2.3.10                
 [55] seqPattern_1.38.0           doParallel_1.0.17          
 [57] codetools_0.2-20            curl_7.0.0                 
 [59] processx_3.8.6              lattice_0.22-7             
 [61] plyr_1.8.9                  Biobase_2.66.0             
 [63] withr_3.0.2                 evaluate_1.0.5             
 [65] Biostrings_2.74.1           pillar_1.11.0              
 [67] MatrixGenerics_1.18.1       whisker_0.4.1              
 [69] KernSmooth_2.23-26          foreach_1.5.2              
 [71] generics_0.1.4              vroom_1.6.5                
 [73] rprojroot_2.1.1             RCurl_1.98-1.17            
 [75] hms_1.1.3                   scales_1.4.0               
 [77] glue_1.8.0                  tools_4.4.2                
 [79] BiocIO_1.16.0               BSgenome_1.74.0            
 [81] GenomicAlignments_1.42.0    fs_1.6.6                   
 [83] XML_3.99-0.18               Cairo_1.6-5                
 [85] impute_1.80.0               plotrix_3.8-4              
 [87] colorspace_2.1-1            GenomeInfoDbData_1.2.13    
 [89] restfulr_0.0.16             cli_3.6.5                  
 [91] S4Arrays_1.6.0              gtable_0.3.6               
 [93] sass_0.4.10                 digest_0.6.37              
 [95] SparseArray_1.6.2           htmlwidgets_1.6.4          
 [97] rjson_0.2.23                farver_2.1.2               
 [99] htmltools_0.5.8.1           lifecycle_1.0.4            
[101] httr_1.4.7                  GlobalOptions_0.1.2        
[103] bit64_4.6.0-1