#####################################################
### rtracklayer: connecting R to genome browsers
#####################################################

### Get islands

## read in aligned reads
library(ShortRead)
setwd("../") # or your path to course data
p <- "extdata/ELAND/080828_HWI-EAS88_0003"
sp <- SolexaPath(p)
filter <- chromosomeFilter("chr2.fa")
aln <- readAligned(sp, "s_1_1_export_head.txt$", 
                   filter = filter)

## calculate coverage and slice into islands
cov <- coverage(aln, start=1)[[1]]
islands <- slice(cov, 8)
islandSums <- viewSums(islands)
max(islandSums)
islands[which.max(islandSums)]

### Convert coverage to a GenomicData

gd <- as(cov, "GenomicData")
names(gd) <- "chr2"

### Accessing a GenomicData

str(start(gd)) ## just the starts
head(as.data.frame(gd)) ## to a data.frame for e.g. plotting

### Subsetting

gd <- gd["chr2"] ## by chromosome
gd <- gd[gd[["score"]] > 10,] ## subsetting

### Convert a GenomicData to an rtracklayer trackSet

library(rtracklayer)
peaksTrack <- trackSet(ranges(gd)[[1]], chrom(gd),
                       dataVals = gd[["score"]],
                       genome = "mm9")

### Output peaks as a WIG file

export(peaksTrack, "peaks.wig", name = "peaks")

### Upload to UCSC manually
## .........

### Upload to UCSC browser automatically

r <- ranges(gd)[[1]]
island <- islands[which.max(islandSums)]
anomaly <- peaksTrack[r %in% island]
session <- browseGenome(anomaly)

## zoom out 4x for context
segment <- genomeSegment(session)
segment <- segment / 4

## display new view
view <- browserView(session, segment)

## display with a cutoff
session <- layTrack(session, anomaly, "anomalyCutoff", yLineMark = 1000,
                    yLineOnOff = TRUE, color = c(0L, 255L, 0L))
view <- browserView(session, full = "anomalyCutoff", hide = "anomaly")

### Get current view coordinates

segment <- genomeSegment(view)

### Download conservation scores in the region

names(tracks(session)) ## list available tracks
track <- trackSet(session, "Conservation",
                  genomeSegment(session) / 100,
                  "phastCons30way")

### GenomicData of conservation scores

GenomicData(IRanges(start(track), end(track)),
            phastCons = dataVals(track), chrom = chrom(segment))

sessionInfo()