###################################################
### chunk number 1: setup
###################################################
options(width=50)
set.seed(0)
library(IRanges)
plotRanges <- function(x, xlim = x, main = deparse(substitute(x)),
                       col = "black", sep = 0.5, ...) 
{
  height <- 1
  if (is(xlim, "Ranges"))
    xlim <- c(min(start(xlim)), max(end(xlim)))
  bins <- disjointBins(IRanges(start(x), end(x) + 1))
  plot.new()
  plot.window(xlim, c(0, max(bins)*(height + sep)))
  ybottom <- bins * (sep + height) - height
  rect(start(x)-0.5, ybottom, end(x)+0.5, ybottom + height, col = col, ...)
  title(main)
  axis(1)
}


###################################################
### chunk number 2: sequences-example
###################################################
s <- c(rep(0, 40), 
       rep(c(1:10, seq(10, 1, -2), rep(1,6), seq(1,5,2), rep(5, 5), 4:1), 
           rep(1:3, 11)),
       rep(0, 50))


###################################################
### chunk number 3: sequences-example-plot
###################################################
plot(s, type="l")


###################################################
### chunk number 4: sequences-diff
###################################################
sum(diff(s) == 0)


###################################################
### chunk number 5: sequences-rle
###################################################
sRle <- Rle(s)
sRle


###################################################
### chunk number 6: rle-basic
###################################################
sRle > 0 | rev(sRle) > 0


###################################################
### chunk number 7: rle-summary
###################################################
sum(sRle > 0)


###################################################
### chunk number 8: rle-cor
###################################################
cor(sRle, rev(sRle))


###################################################
### chunk number 9: rle-split
###################################################
sRleList <- split(sRle, rep(c("chr1", "chr2"), each = length(sRle)/2))
sRleList[1]


###################################################
### chunk number 10: ranges-construct
###################################################
ir <- IRanges(c(1, 8, 14, 15, 19, 34, 40), 
              width=c(12, 6,  6, 15, 6, 2, 7))


###################################################
### chunk number 11: ranges-plot
###################################################
plotRanges(ir)


###################################################
### chunk number 12: ranges-start
###################################################
start(ir)


###################################################
### chunk number 13: ranges-end
###################################################
end(ir)


###################################################
### chunk number 14: ranges-width
###################################################
width(ir)


###################################################
### chunk number 15: ranges-subset
###################################################
ir[1:5]


###################################################
### chunk number 16: ranges-split
###################################################
rl <- split(ir, c(rep("chr1", 2), rep("chr2", 3), "chr1", "chr2"))
rl[1]


###################################################
### chunk number 17: ranges-shift
###################################################
shift(ir, 1)


###################################################
### chunk number 18: ranges-shift-plot
###################################################
plotRanges(shift(ir, 1))


###################################################
### chunk number 19: ranges-resize
###################################################
ir15 <- resize(ir, 15) # forward
ir15 <- resize(ir, 15, start=FALSE) # backward
ir15[1:3]


###################################################
### chunk number 20: ranges-restrict
###################################################
restrict(ir15, 1)


###################################################
### chunk number 21: ranges-reduce
###################################################
reduce(ir)


###################################################
### chunk number 22: ranges-reduce-plot
###################################################
plotRanges(reduce(ir))


###################################################
### chunk number 23: ranges-gaps
###################################################
gaps(ir)


###################################################
### chunk number 24: ranges-gaps-plot
###################################################
plotRanges(gaps(ir), ir)


###################################################
### chunk number 25: ranges-disjoin
###################################################
disjoin(ir)


###################################################
### chunk number 26: ranges-disjoin-plot
###################################################
plotRanges(disjoin(ir))


###################################################
### chunk number 27: ranges-overlap
###################################################
ol <- findOverlaps(ir, reduce(ir))
as.matrix(ol)


###################################################
### chunk number 28: ranges-coverage
###################################################
cov <- coverage(ir)


###################################################
### chunk number 29: ranges-coverage-plot
###################################################
plotRanges(ir)
cov <- as.vector(cov)
mat <- cbind(seq_along(cov)-0.5, cov)
d <- diff(cov) != 0
mat <- rbind(cbind(mat[d,1]+1, mat[d,2]), mat)
mat <- mat[order(mat[,1]),]
lines(mat, col="red", lwd=4)
axis(2)


###################################################
### chunk number 30: ranges-coverage
###################################################
covL <- coverage(rl)
covL


###################################################
### chunk number 31: views-sequence-plot
###################################################
plot(s, type="l")
abline(h = 3, col="red")


###################################################
### chunk number 32: views-rle
###################################################
Views(sRle, as(sRle > 3, "IRanges"))


###################################################
### chunk number 33: views-slice
###################################################
sViews <- slice(sRle, 4)


###################################################
### chunk number 34: views-slice-list
###################################################
sViewsList <- slice(sRleList, 4)
sViewsList[1]


###################################################
### chunk number 35: views-viewSums
###################################################
viewSums(sViews)
viewSums(sViewsList)


###################################################
### chunk number 36: views-viewMaxs
###################################################
viewMaxs(sViews)
viewMaxs(sViewsList)


###################################################
### chunk number 37: naive-RangedData-1
###################################################
chr <- c("chr1", "chr2", "chr1")
strand <- c("+", "+", "-")
start <- c(3L, 4L, 1L)
end <- c(7L, 5L, 3L)
score <- c(1L, 3L, 2L)
naiveTable <- data.frame(chr, strand, score, start, end)
naiveTable


###################################################
### chunk number 38: naive-RangedData-2
###################################################
getRange <- function(x) range(x[c("start", "end")])
by(naiveTable, naiveTable[["chr"]], getRange)


###################################################
### chunk number 39: RangedData-construction
###################################################
rdTable <- RangedData(IRanges(start, end), strand, score, space = chr)


###################################################
### chunk number 40: RangedData-display
###################################################
rdTable


###################################################
### chunk number 41: rd-ranges
###################################################
ranges(rdTable)[1]


###################################################
### chunk number 42: rd-values
###################################################
values(rdTable)[1]


###################################################
### chunk number 43: attribute-accessors
###################################################
start(rdTable)


###################################################
### chunk number 44: access-column
###################################################
rdTable$strand


###################################################
### chunk number 45: subset-features
###################################################
## get the first 2 features
first10 <- rdTable[1:2,]
## get features on positive strand
pos <- rdTable[rdTable$strand == "+",]
## get first chromosome
chr1Table <- rdTable[1]
## with only the score column
scoreTable <- rdTable[,"score"]


###################################################
### chunk number 46: rle-to-rd
###################################################
head(as(sRleList, "RangedData"), 3)


###################################################
### chunk number 47: rle-views-rd
###################################################
height <- unlist(viewMaxs(sViewsList))
RangedData(sViewsList, height)


###################################################
### chunk number 48: rd-to-rle
###################################################
coverage(rdTable, weight = "score")[1]