In this analysis we are going to calculate correlations of the PARP1 reads with a variety of histone mark ChIP-seq data from UCSC/ENCODE. The PARP1 reads and histone mark data has already been processed and is made available as part of the fmdatabreastcaparp1 package.
knitr::opts_chunk$set(dev='png')
library(GenomicRanges) library(ggplot2) library(fmcorrelationbreastcaparp1) library(fmdatabreastcaparp1) library(BSgenome.Hsapiens.UCSC.hg19)
data(parp1_ln4_unique) data(parp1_ln5_unique) data(histone_marks) data(tss_windows)
Calculate the weighted coverage of the PARP1 mcf7 and mdamb231 sample reads, and then sum the reads in each TSS window.
mcf7_cov <- coverage(parp1_ln4_unique, weight = "n_count") mdamb231_cov <- coverage(parp1_ln5_unique, weight = "n_count") tss_windows <- binned_function(tss_windows, mcf7_cov, "sum", "parp1_mcf7") tss_windows <- binned_function(tss_windows, mdamb231_cov, "sum", "parp1_mdamb231")
Get the averaged ChIP-seq peak intensity in each tss window.
for (i_name in names(histone_marks)){ histone_cov <- coverage(histone_marks[[i_name]], weight = "mcols.signal") tss_windows <- binned_function(tss_windows, histone_cov, "mean_nozero", i_name) }
Now with the Parp1 reads and histone mark signal added to the TSS's, we can start doing some correlations.
non_zero <- "both"
h3k4me3k_v_mcf7 <- subsample_nonzeros(mcols(tss_windows), c("H3k4me3_r1", "parp1_mcf7"), non_zero = non_zero, n_points = 10000) ggplot(h3k4me3k_v_mcf7, aes(x = H3k4me3_r1, y = parp1_mcf7)) + geom_point() + scale_y_log10() + scale_x_log10() cor(log10(h3k4me3k_v_mcf7[,1]+1), log10(h3k4me3k_v_mcf7[,2]+1)) h3k27ac_v_mcf7 <- subsample_nonzeros(mcols(tss_windows), c("H3k27ac", "parp1_mcf7"), non_zero = non_zero, n_points = 10000) ggplot(h3k27ac_v_mcf7, aes(x = H3k27ac, y = parp1_mcf7)) + geom_point() + scale_y_log10() + scale_x_log10() cor(log10(h3k27ac_v_mcf7[,1]+1), log10(h3k27ac_v_mcf7[,2]+1))
Cool. Now we are showing some promise. Let's do them all.
all_comb <- expand.grid(names(histone_marks), c("parp1_mcf7", "parp1_mdamb231"), stringsAsFactors = FALSE) out_cor <- lapply(seq(1, nrow(all_comb)), function(i_row){ #print(i_row) correlate_non_zero(mcols(tss_windows), as.character(all_comb[i_row,]), log_transform = TRUE, non_zero = non_zero, test = TRUE) }) all_comb_names <- paste(all_comb[,1], all_comb[,2], sep = "_v_") out_cor <- do.call(rbind, out_cor) rownames(out_cor) <- all_comb_names
TSS correlations:
knitr::kable(out_cor)
out_graphs <- lapply(seq(1, nrow(all_comb)), function(i_row){ use_vars <- as.character(all_comb[i_row,]) subpoints <- subsample_nonzeros(mcols(tss_windows), use_vars, non_zero = non_zero, n_points = 10000) ggplot(subpoints, aes_string(x = use_vars[1], y = use_vars[2])) + geom_point() + scale_y_log10() + scale_x_log10() }) out_graphs
Are these correlations a result of association with the TSS's? One way to test this is to set up a calculation genome-wide.
genome_tiles <- tileGenome(seqinfo(Hsapiens), tilewidth = 2000, cut.last.tile.in.chrom = TRUE) genome_tiles <- binned_function(genome_tiles, mcf7_cov, "sum", "parp1_mcf7") genome_tiles <- binned_function(genome_tiles, mdamb231_cov, "sum", "parp1_mdamb231")
for (i_name in names(histone_marks)){ histone_cov <- coverage(histone_marks[[i_name]], weight = "mcols.signal") genome_tiles <- binned_function(genome_tiles, histone_cov, "mean_nozero", i_name) }
genome_h3k4me3k_v_mcf7 <- subsample_nonzeros(mcols(genome_tiles), c("H3k4me3_r1", "parp1_mcf7"), non_zero = non_zero, n_points = 10000) ggplot(genome_h3k4me3k_v_mcf7, aes(x = H3k4me3_r1, y = parp1_mcf7)) + scale_x_log10() + scale_y_log10() + geom_point() cor(log(genome_h3k4me3k_v_mcf7[,1]+1), log(genome_h3k4me3k_v_mcf7[,2]+1))
genome_cor <- lapply(seq(1, nrow(all_comb)), function(i_row){ #print(i_row) correlate_non_zero(mcols(genome_tiles), as.character(all_comb[i_row,]), log_transform = TRUE, non_zero = non_zero, test = TRUE) }) all_comb_names <- paste(all_comb[,1], all_comb[,2], sep = "_v_") genome_cor <- do.call(rbind, genome_cor) rownames(genome_cor) <- all_comb_names
Genome wide correlations:
knitr::kable(genome_cor)
We will save the correlation results in some plain text files.
saveloc <- "../inst/correlation_tables" write.table(out_cor, file = file.path(saveloc, "histone_marks_tss.txt"), sep = "\t") write.table(genome_cor, file = file.path(saveloc, "histone_marks_genome.txt"), sep = "\t")
Sys.time() sessionInfo()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.