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