################################################### ### chunk number 1: eval=FALSE ################################################### ## library( ShortRead ) ## maps.me1 <- sapply( list.files( "ShortReadExampleData/H3K4me1", "run.*lane.\\.map" ), ## function(filename) ## readAligned( "ShortReadExampleData/H3K4me1", filename, type="MAQMapShort" ) ) ## maps.me3 <- sapply( list.files( "ShortReadExampleData/H3K4me3", "run.*lane.\\.map" ), ## function(filename) ## readAligned( "ShortReadExampleData/H3K4me3", filename, type="MAQMapShort" ) ) ## ## get_coverage <- function( lanes, chrom, chromLength ) { ## res <- Rle( 0, chromLength ) ## for ( lane in lanes ) { ## lc <- lane[ chromosome(lane) == chrom ] ## cvg <- coverage( lc, start = 1, end = chromLength, ## extend = 185L - width(lc) ) ## res <- res + cvg[[ chrom ]] ## } ## res ## } ## ## lengthChr10 <- 135374737 ## ## coverage.chr10.me1 <- get_coverage( maps.me1, "10", lengthChr10 ) ## coverage.chr10.me3 <- get_coverage( maps.me3, "10", lengthChr10 ) ## ## save( coverage.chr10.me1, coverage.chr10.me3, file="coverage.H3K4meX.rda") ################################################### ### chunk number 2: load ################################################### library("ShortRead") load(file.path("..", "data", "coverage.H3K4meX.rda")) ################################################### ### chunk number 3: coverage ################################################### coverage.chr10.me1 coverage.chr10.me3 ################################################### ### chunk number 4: toRangedData ################################################### library("rtracklayer") coverage.chr10.me1.rd <- as( coverage.chr10.me1, "RangedData" ) coverage.chr10.me3.rd <- as( coverage.chr10.me3, "RangedData" ) names( coverage.chr10.me1.rd ) <- "chr10" names( coverage.chr10.me3.rd ) <- "chr10" ################################################### ### chunk number 5: rtracklayer ################################################### export( coverage.chr10.me1.rd, "chr10_me1.wig", format="wig" ) export( coverage.chr10.me3.rd, "chr10_me3.wig", format="wig" ) ################################################### ### chunk number 6: rtracklayerWorkAround ################################################### names( coverage.chr10.me1.rd ) <- "10" names( coverage.chr10.me3.rd ) <- "10" export( coverage.chr10.me1.rd, "chr10_me1.wig", format="wig" ) export( coverage.chr10.me3.rd, "chr10_me3.wig", format="wig" ) ################################################### ### chunk number 7: hilbert eval=FALSE ################################################### ## library( HilbertVisGUI ) ## hilbertDisplay( as.integer( coverage.chr10.me1 ) ) ################################################### ### chunk number 8: hilbert2 eval=FALSE ################################################### ## hilbertDisplay( ## as.integer( coverage.chr10.me1 ), ## as.integer( coverage.chr10.me3 ) ) ################################################### ### chunk number 9: subseq ################################################### as.vector( subseq( coverage.chr10.me1, 123456-50, 123456+50 ) ) ################################################### ### chunk number 10: deltaConvolve ################################################### deltaConvolve <- function( coverage, winCentres, revStrand, left, right) { result <- rep( 0, right - left + 1 ) for( i in seq(along=winCentres) ) { v <- as.vector( subseq( coverage, winCentres[i] + left, winCentres[i] + right ) ) if( revStrand[i] ) v <- rev( v ) result <- result + v } result } ################################################### ### chunk number 11: biomaRt ################################################### library("biomaRt") mart <- useMart( "ensembl", dataset = "hsapiens_gene_ensembl" ) martReply <- getBM( attributes=c( "transcript_start", "transcript_end", "strand" ), filters="chromosome_name", values="10", mart=mart ) ################################################### ### chunk number 12: martReply ################################################### head(martReply) ################################################### ### chunk number 13: tss ################################################### tss <- ifelse( martReply$strand == 1, martReply$transcript_start, martReply$transcript_end ) ################################################### ### chunk number 14: calldc ################################################### win = c(-500, 500) H3K4me1 <- deltaConvolve( coverage.chr10.me1, tss, martReply$strand == -1, left=win[1], right=win[2] ) H3K4me3 <- deltaConvolve( coverage.chr10.me3, tss, martReply$strand == -1, left=win[1], right=win[2] ) ################################################### ### chunk number 15: plotfun ################################################### plotfun = function( win, ... ) { y = cbind( ... ) matplot( win[1]:win[2], y, type='l', lty=1, lwd=2, ylab="average coverage", xlab="distance from TSS" ) legend( "topleft", colnames(y), fill=seq_len(ncol(y)) ) } ################################################### ### chunk number 16: plot1 ################################################### plotfun( win, H3K4me1, H3K4me3 ) ################################################### ### chunk number 17: normalise ################################################### H3K4me1.N = H3K4me1 / sum(coverage.chr10.me1) H3K4me3.N = H3K4me3 / sum(coverage.chr10.me3) ################################################### ### chunk number 18: plot2 ################################################### plotfun( win, H3K4me1.N, H3K4me3.N ) ################################################### ### chunk number 19: wideWindow ################################################### win.w = c(-3000, 3000) H3K4me1.w <- deltaConvolve(coverage.chr10.me1, tss, martReply$strand == -1, left=win.w[1], right=win.w[2]) H3K4me3.w <- deltaConvolve(coverage.chr10.me3, tss, martReply$strand == -1, left=win.w[1], right=win.w[2]) H3K4me1.w.N = H3K4me1.w / sum(coverage.chr10.me1) H3K4me3.w.N = H3K4me3.w / sum(coverage.chr10.me3) plotfun( win.w, H3K4me1.w.N, H3K4me3.w.N )