### R code from vignette source 'LiebermanAidenHiC2009.Rnw' ################################################### ### code chunk number 1: data ################################################### library(LiebermanAidenHiC2009) data("HiC_GM_chr14") head(HiC_GM_chr14) pos = with(HiC_GM_chr14, cbind(position1, position2)) ################################################### ### code chunk number 2: figscatter1 (eval = FALSE) ################################################### ## plot(pos, pch='.', col="#77777777") ################################################### ### code chunk number 3: figscatter2 ################################################### png("figscatter.png",width=1024,height=1024) plot(pos, pch='.', col="#77777777") dev.off() ################################################### ### code chunk number 4: smooth ################################################### chrlen = max(pos) gridsize = ceiling(chrlen/2e5) bandwidth = 3e5 den = bkde2D( pos, bandwidth=c(1,1)*bandwidth, gridsize=c(1,1)*gridsize) ################################################### ### code chunk number 5: symm ################################################### den$fhat <- den$fhat + t(den$fhat) ################################################### ### code chunk number 6: figmatrix1 ################################################### with(den, image(x=x1, y=x2, z=fhat^0.3, col=colorRampPalette(c("white","blue"))(256), useRaster=TRUE)) ################################################### ### code chunk number 7: normalize ################################################### m = matrix(0, nrow=gridsize, ncol=gridsize) for(i in 1:(gridsize-1)) { band = (row(m)==col(m)+i) m[band] = mean(den$fhat[band]) } m = m + t(m) diag(m) = mean(diag(den$fhat)) ################################################### ### code chunk number 8: figmatm ################################################### image(x=den$x1, y=den$x2, z=m^0.3, col=colorRampPalette(c("white","blue"))(256), useRaster=TRUE) ################################################### ### code chunk number 9: figaverage ################################################### genomicDistance = den$x1 - min(den$x1) averageInteractions = m[1,] plot(genomicDistance, sqrt(averageInteractions), type ="l") ################################################### ### code chunk number 10: normalize ################################################### fhatNorm <- den$fhat/m ################################################### ### code chunk number 11: figmatrix2 ################################################### image(x=den$x1, y=den$x2, z=fhatNorm, col=colorRampPalette(c("white","blue"))(256), useRaster=TRUE) ################################################### ### code chunk number 12: cormatrix ################################################### cm <- cor(fhatNorm) ################################################### ### code chunk number 13: figmatrix3 ################################################### image(x=den$x1, y=den$x2, z=cm, zlim=c(-1,1), col=colorRampPalette(c("red", "white","blue"))(256), useRaster=TRUE) ################################################### ### code chunk number 14: cormatrix ################################################### princp = princomp(cm) ################################################### ### code chunk number 15: figprcomp1 ################################################### plot(den$x1, princp$loadings[,1], type="l") ################################################### ### code chunk number 16: pc1Vec ################################################### pc1Vec = Rle(values = princp$loadings[,1], lengths = c(den$x1[1], diff(den$x1))) pc2Vec = Rle(values = princp$loadings[,2], lengths = c(den$x1[1], diff(den$x1))) ################################################### ### code chunk number 17: readChipseq ################################################### createRleVector = function(tab){ RleVec = Rle(0, max(tab$end)) for(i in 1:nrow(tab)){ RleVec = RleVec + Rle(values = c(0, tab$signalValue[i], 0), lengths = c(tab$start[i]-1, tab$end[i]-tab$start[i]+1, length(RleVec)-tab$end[i])) } RleVec } ################################################### ### code chunk number 18: read ################################################### data("ChipSeqData") H3K27me3 = createRleVector(H3K27me3.df) H3K36me3 = createRleVector(H3K36me3.df) DNAse1 = createRleVector(DNAse1.df) DNAse2 = createRleVector(DNAse2.df) ################################################### ### code chunk number 19: combinednase ################################################### length(DNAse1) length(DNAse2) DNAse = DNAse1 + DNAse2[seq(along=DNAse1)] ################################################### ### code chunk number 20: plotRle ################################################### plotRle = function(RleVector, ...){ plot(end(RleVector), runValue(RleVector)+1, type="h", log="y", xlim = c(1.5e+7, 107000000), xlab="", ylab=deparse(substitute(RleVector)), ...) } ################################################### ### code chunk number 21: figrle ################################################### par(mfrow=c(4,1), mai=c(0.5,0.7,0.1,0.1)) plotRle(pc1Vec) plotRle(H3K27me3) plotRle(H3K36me3) plotRle(DNAse1) ################################################### ### code chunk number 22: correlation ################################################### c(length(H3K27me3), length(H3K36me3), length(DNAse), length(pc1Vec)) x = seq(along=H3K36me3) cor(H3K27me3[x], pc1Vec[x]) cor(H3K36me3, pc1Vec[x]) cor(DNAse[x], pc1Vec[x]) ################################################### ### code chunk number 23: pData (eval = FALSE) ################################################### ## sampleAnnotation = ## data.frame( ## file = c("GSM455133_30E0LAAXX.1.maq.hic.summary.binned.txt.gz", ## "GSM455134_30E0LAAXX.2.maq.hic.summary.binned.txt.gz", ## "GSM455135_30U85AAXX.2.maq.hic.summary.binned.txt.gz", ## "GSM455136_30U85AAXX.3.maq.hic.summary.binned.txt.gz", ## "GSM455137_30305AAXX.1.maq.hic.summary.binned.txt.gz", ## "GSM455138_30305AAXX.2.maq.hic.summary.binned.txt.gz"), ## experiment = c(1, 1, 2, 2, 3, 3), ## lane = c(1,2,1,2,1,2), ## restrictionenzyme = c("HindIII", "HindIII", "HindIII", "HindIII", ## "NcoI", "NcoI"), ## stringsAsFactors = FALSE) ################################################### ### code chunk number 24: readtable (eval = FALSE) ################################################### ## inputDir = "WHERE YOU SAVED THE FILES" ## outputDir = "WHERE YOU WANT TO SAVE THE FILES" ## df = vector(mode="list", length=nrow(sampleAnnotation)) ## for(i in seq(along=df)) { ## cat("Reading file", i, "\n") ## r = read.table(gzfile(file.path(inputDir, sampleAnnotation$file[i])), ## header=FALSE, sep="\t", comment.char = "", stringsAsFactors=FALSE) ## colnames(r) = c("read name", ## "chromosome1", "position1", "strand1", "restrictionfragment1", ## "chromosome2", "position2", "strand2", "restrictionfragment2") ## ## filter chromosome 14 ## df[[i]] = subset(r, (chromosome1==14L) & (chromosome2==14L)) ## } ## HiC_GM_chr14 = do.call(rbind, df) ## save(HiC_GM_chr14, ## file=file.path(outputDir, "HiC_GM_chr14.RData")) ################################################### ### code chunk number 25: ChIPfiles (eval = FALSE) ################################################### ## files = c("wgEncodeBroadChipSeqPeaksGm12878H3k27me3.broadPeak.gz", ## "wgEncodeBroadChipSeqPeaksGm12878H3k36me3.broadPeak.gz", ## "wgEncodeUwDnaseSeqHotspotsRep1Gm06990.broadPeak.gz", ## "wgEncodeUwDnaseSeqHotspotsRep2Gm06990.broadPeak.gz") ################################################### ### code chunk number 26: readtable (eval = FALSE) ################################################### ## inputDir = "WHERE YOU SAVED THE FILES" ## outputDir = "WHERE YOU WANT TO SAVE THE FILES" ## cs = vector(mode="list", length=length(files)) ## for(i in seq(along=files)) { ## cat("Reading file", i, "\n") ## tab = read.table(gzfile(file.path(inputDir, files[i])), header=FALSE, sep="\t", comment.char = "", stringsAsFactors=FALSE) ## colnames(tab) <- c("chr", "start", "end", "name", "score", "strand", ## "signalValue", "pValue", "qValue") ## cs[[i]] = subset(tab, chr=="chr14") ## } ## ## H3K27me3.df = cs[[1]] ## H3K36me3.df = cs[[2]] ## DNAse1.df = cs[[3]] ## DNAse2.df = cs[[4]] ## ## save(list=c("H3K27me3.df", "H3K36me3.df", "DNAse1.df", "DNAse2.df"), file=file.path(outputDir, "ChipSeqData.RData")) ################################################### ### code chunk number 27: sessionInfo ################################################### toLatex(sessionInfo())