library(fastcluster) library(DESeq) source("../script/heatmap_2_multibar.R") color.ramp.size <- 10 load("vsd.RData") clust.cols <- hclust(as.dist(read.table("vsd_tf_cols.phylip", skip=1, fill=T, row.names=1, col.names=c(NA, 1:scan("vsd_tf_cols.phylip", nlines=1)))), method="average") clust.rows <- hclust(as.dist(read.table("vsd_tf_rows.phylip.gz", skip=1, fill=T, row.names=1, col.names=c(NA, 1:scan("vsd_tf_rows.phylip.gz", nlines=1)))), method="average") attach(envs[["tf"]]) cell.type <- factor(sapply(strsplit(colnames(vsd), " "), function(x) x[1]), levels=rev(c("HeLa-S3", "H1-hESC", "HepG2", "K562", "GM12878"))) antibody <- factor(sapply(strsplit(colnames(vsd), " "), 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 <- estimateDispersions(estimateSizeFactors(newCountDataSet(counts.pol2, expt.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.color.ramp <- colorRampPalette(c("white", "red"))(color.ramp.size) 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) } }) # normalize histone data 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) cells <- rev(c("HeLa-S3", "H1-hESC", "HepG2", "K562", "GM12878")) color.table <- rainbow(nlevels(cell.type)) names(color.table) <- levels(cell.type) save.image("../figure/all_cell_heatmap.RData") png("../figure/allCellHeatmap.png", width=4, height=6, units="in", res=300) 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), labRow=NA, labCol=NA, trace="none", cexCol=2, key="F", margins=c(0,0), lwid=c(1E-10,1.5,1), ColSideColors=color.table[cell.type], RowSideColors=color.mat, lwd=0.01) dev.off()