###################################################
### chunk number 1: setup
###################################################
options(width = 40)


###################################################
### chunk number 2: peaks
###################################################
library(ShortRead)
load("../data/ctcf.rda")
library(BSgenome.Mmusculus.UCSC.mm9)
mouse.chromlens <- seqlengths(Mmusculus)
cov.ctcf <- coverage(ctcf, width = mouse.chromlens, extend = 126L)
peaks <- slice(cov.ctcf$chr10, lower = 8)
peak.depths <- viewMaxs(peaks)
peak.areas <- viewSums(peaks)
wpeaks <- tail(order(peak.depths), 4)
library(rtracklayer)


###################################################
### chunk number 3: create-track
###################################################
peakTrack <- GenomicData(peaks, depth = peak.depths, 
                         area = peak.areas, 
                         chrom = "chr10", 
                         genome = "hg18")
peakTrack


###################################################
### chunk number 4: attribute-accessors
###################################################
head(start(peakTrack))


###################################################
### chunk number 5: sol1-strand-accessor
###################################################
head(width(peakTrack))


###################################################
### chunk number 6: sol2-width-accessor
###################################################
genome(peakTrack)


###################################################
### chunk number 7: access-column
###################################################
head(peakTrack$area)


###################################################
### chunk number 8: sol3-reconstruct
###################################################
data.frame(start = start(peakTrack), 
           width = width(peakTrack),
           area = peakTrack$area,
           depth = peakTrack$depth)
## easier:
as.data.frame(peakTrack)


###################################################
### chunk number 9: subset-features
###################################################
## get the first 10 targets
first10 <- peakTrack[1:10,]
## get peaks of depth > 12
highPeaks <- peakTrack[peakTrack$depth > 12,]
## get first (and only) chromosome
chrPeaks <- peakTrack[1]


###################################################
### chunk number 10: sol4-subset
###################################################
peakTrack[peakTrack$area > 1000,"area"]


###################################################
### chunk number 11: export
###################################################
export(peakTrack, "peaks.bed")
restoredTrack <- import("peaks.bed")
## as character vector
peaksChar <- export(peakTrack, format = "gff1")


###################################################
### chunk number 12: sol5-export-gff
###################################################
export(peakTrack, "peaks.gff")


###################################################
### chunk number 13: sol6-import-gff
###################################################
peakGff <- import("peaks.gff", 
                    genome = "hg18")


###################################################
### chunk number 14: browserSession
###################################################
session <- browserSession("UCSC")


###################################################
### chunk number 15: genomeBrowsers
###################################################
genomeBrowsers()


###################################################
### chunk number 16: layTrack
###################################################
track(session, "peaks") <- peakTrack
## equivalently
session$peaks <- peakTrack


###################################################
### chunk number 17: sol7-load-track eval=FALSE
###################################################
## session$topPeaks <- peakTrack[wpeaks,]


###################################################
### chunk number 18: first-site-range
###################################################
high <- tail(wpeaks, 1)
region <- range(peakTrack[high,])


###################################################
### chunk number 19: first-site-zoom
###################################################
region <- region * -10


###################################################
### chunk number 20: coverageTrack
###################################################
covTrack <- as(cov.ctcf$chr10, "RangedData")
names(covTrack) <- "chr10"
session$coverage <- covTrack


###################################################
### chunk number 21: browserView
###################################################
view <- browserView(session, region, 
                    full = "coverage")


###################################################
### chunk number 22: sleep
###################################################
Sys.sleep(20)


###################################################
### chunk number 23: sol8-create-view
###################################################
viewOut <- browserView(session, range(view) * -2)


###################################################
### chunk number 24: browseGenome
###################################################
browseGenome(peakTrack, 
             range = range(peakTrack[high,]) * -10)


###################################################
### chunk number 25: set-view-range
###################################################
## zoom in 2X
range(view) <- range(view) * 2 


###################################################
### chunk number 26: sol9-view-second-site
###################################################
second <- tail(wpeaks, 2)[1]
range(viewOut) <- range(peakTrack[second,]) * -5


###################################################
### chunk number 27: visible
###################################################
## hide the Conservation track
visible(view)["Conservation"] <- FALSE


###################################################
### chunk number 28: sol10-show-conservation
###################################################
visible(view)["Ensembl Genes"] <- TRUE


###################################################
### chunk number 29: get-track-names
###################################################
head(trackNames(session))


###################################################
### chunk number 30: download-track
###################################################
cons <- track(session, "Conservation")
## or specific region
cons <- track(session, "Conservation", 
              range(view) * 2)
## shortcut
session$Conservation


###################################################
### chunk number 31: session-info
###################################################
sessionInfo()