###################################################
### 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 )