library(fastcluster) library(DESeq) source("../script/heatmap_2_multibar.R") color.ramp.size <- 10 load("vsd.RData") cells <- rev(c("HeLa-S3", "H1-hESC", "HepG2", "K562", "GM12878")) color.table <- rainbow(length(cells)) names(color.table) <- cells pol2.color.ramp <- colorRampPalette(c("white", "red"))(color.ramp.size) for (cell in cells) { with(envs[[cell]], { clust.cols <- hclust(as.dist(read.table(paste("vsd_", cell, "_cols.phylip", sep=""), skip=1, fill=T, row.names=1, col.names=c(NA, 1:scan(paste("vsd_", cell, "_cols.phylip", sep=""), nlines=1)))), method="average") clust.rows <- hclust(as.dist(read.table(paste("vsd_", cell, "_rows.phylip.gz", sep=""), skip=1, fill=T, row.names=1, col.names=c(NA, 1:scan(paste("vsd_", cell, "_rows.phylip.gz", sep=""), nlines=1)))), method="average") transfac <- sapply(strsplit(colnames(counts), " "), function(x) x[2]) # normalize Pol2 data counts.pol2 <- regions[regions[,"nearest Pol2 peak"] != '', grep("pol2.+rep", colnames(regions))] expt.pol2 <- factor(sub(" rep.*", "", colnames(counts.pol2))) cds.pol2 <- newCountDataSet(counts.pol2, expt.pol2) cds.pol2 <- estimateSizeFactors(cds.pol2) cds.pol2 <- estimateDispersions(cds.pol2, method="pooled", sharingMode="fit-only", fitType="local") # pooled so they can be meaned vsd.pol2 <- getVarianceStabilizedData(cds.pol2) pol2 <- rowMeans(vsd.pol2) pol2.min <- min(pol2) pol2.range <- max(pol2) - pol2.min # make vector of colors pol2.colors <- sapply(rownames(regions), function(region) { if (regions[region,"nearest Pol2 peak"] == "") pol2.color.ramp[1] else { # NA -> 0 pol2.color.ramp[round((length(pol2.color.ramp) - 1) * (pol2[region] - pol2.min) / pol2.range) + 1] # normalize to the range 1:length(pol2.color.ramp) to pick element of color ramp (this is not the elegant way to do it) } }) counts.h3k4 <- regions[,grep("H3K4", colnames(regions))] expt.h3k4 <- factor(sub(" Bernstein.*", "", colnames(counts.h3k4))) cds.h3k4 <- estimateDispersions(estimateSizeFactors(newCountDataSet(counts.h3k4, expt.h3k4)), method="pooled", sharingMode="fit-only", fitType="local") # pooled so they can be meaned vsd.h3k4 <- getVarianceStabilizedData(cds.h3k4) # make matrix of colors me1 <- rowMeans(vsd.h3k4[,grep("H3K4me1", colnames(vsd.h3k4))]) me1.color.ramp <- colorRampPalette(c("white", "purple"))(color.ramp.size) me1.min <- min(me1) me1.range <- max(me1) - me1.min me1.colors <- me1.color.ramp[round((length(me1.color.ramp) - 1) * (me1 - me1.min) / me1.range) + 1] me3 <- rowMeans(vsd.h3k4[,grep("H3K4me3", colnames(vsd.h3k4))]) me3.color.ramp <- colorRampPalette(c("white", "darkgreen"))(color.ramp.size) me3.min <- min(me3) me3.range <- max(me3) - me3.min me3.colors <- me3.color.ramp[round((length(me3.color.ramp) - 1) * (me3 - me3.min) / me3.range) + 1] color.mat <- cbind(me1.colors, me3.colors, pol2.colors) }) } totalcols <- sum(sapply(cells, function(cell) with(envs[[cell]], ncol(vsd)))) save.image("cell_specific_heatmaps.RData") for (cell in cells) { png(paste0("cellSpecificHeatmap", cell, ".png"), height=11-2*1, width=8.5-2*1, units="in", res=300) with(envs[[cell]], { heatmap.2.multibar(vsd, Rowv=as.dendrogram(clust.rows), Colv=as.dendrogram(clust.cols), dendrogram="c", scale="none", col=colorRampPalette(c("white", "darkblue"))(color.ramp.size), labCol=transfac, labRow=NA, trace="none", key=F, margins=c(5,0), lwid=c(1E-10,3,1), RowSideColors=color.mat, main="", lwd=0.01) }) dev.off() }