### R code from vignette source 'vignettes/ggbio/inst/doc/intro.Rnw'

###################################################
### code chunk number 1: options
###################################################
options(width=72)
## theme_set(theme_grey(base_size = 9))


###################################################
### code chunk number 2: loading
###################################################
library(ggbio)


###################################################
### code chunk number 3: ggbio-loading
###################################################
p <- qplot(data = mtcars, mpg, cyl)
print(p)


###################################################
### code chunk number 4: GRanges-sample
###################################################
set.seed(1)
N <- 1000
library(GenomicRanges)
gr <- GRanges(seqnames = 
              sample(c("chr1", "chr2", "chr3"),
                     size = N, replace = TRUE),
              IRanges(
                      start = sample(1:300, size = N, replace = TRUE),
                      width = sample(70:75, size = N,replace = TRUE)),
              strand = sample(c("+", "-", "*"), size = N, 
                replace = TRUE),
              value = rnorm(N, 10, 3), score = rnorm(N, 100, 30),
              group = sample(c("Normal", "Tumor"), 
                size = N, replace = TRUE),
              pair = sample(letters, size = N, 
                replace = TRUE))


###################################################
### code chunk number 5: qplot-gr-full
###################################################
p <- qplot(gr)
print(p)


###################################################
### code chunk number 6: qplot-gr-full-themebw
###################################################
p <- p + theme_bw()
print(p)


###################################################
### code chunk number 7: qplot-gr-cov-line
###################################################
p <- qplot(gr, nrow = 1, geom = "coverage.p")
p <- p + geom_hline(yintercept = 40, color = "red", 
    size = 1)
print(p)


###################################################
### code chunk number 8: qplot-gr-all
###################################################
gr.sub <- gr[seqnames(gr) == "chr1"] #or 
p1 <- qplot(gr.sub, geom = "full") + opts(title = "full")
p2 <- qplot(gr.sub, geom = "point", y = value) + opts(title = "point")
p3 <- qplot(gr.sub, geom = "line", y = value) + opts(title = "line")
p4 <- qplot(gr.sub, geom = "coverage.line") + opts(title = "coverage.line")   
p5 <- qplot(gr.sub, geom = "coverage.polygon") + opts(title = "coverage.polygon")   
library(gridExtra)
grid.arrange(p1, p2, p3, p4, p5,ncol = 2)


###################################################
### code chunk number 9: gr-ncol2
###################################################
p <- qplot(gr, ncol = 2)
print(p)


###################################################
### code chunk number 10: facet-group
###################################################
  p <- qplot(gr, facets = group ~ seqnames)
  print(p)


###################################################
### code chunk number 11: facet-gr
###################################################
 gr.region <- GRanges(c("chr1", "chr2", "chr3"), 
  IRanges(c(100, 200, 250), 
  width = 70))
## facet_grid
p <- qplot(gr, facet_gr = gr.region)
## facet_wrap
p <- qplot(gr, facet_gr = gr.region, nrow = 2) + 
  scale_y_continuous(limits = c(0, 90))
print(p)


###################################################
### code chunk number 12: chevron1
###################################################
gr.chr1 <- GRanges("chr1", IRanges(c(100, 200, 300), width = 50))
p <- qplot(gr.chr1)
gr.gaps <- gaps(gr.chr1)[-1]
values(gr.gaps)$score <- c(1, 100)
p1 <- p + geom_chevron(gr.gaps)
p2 <- p + geom_chevron(gr.gaps, aes(size = score), offset = "score",
                 chevron.height = c(0.1, 0.2))
p3 <- p + geom_chevron(gr.gaps, offset = -0.1)
tracks(p1, p2, p3)


###################################################
### code chunk number 13: splice-grl (eval = FALSE)
###################################################
## library(TxDb.Hsapiens.UCSC.hg19.knownGene)
## data(genesymbol)
## txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
## exons.tx <- exonsBy(txdb, by = "tx")
## exons.rbm17 <- subsetByOverlaps(exons.tx, genesymbol["RBM17"])
## nms <- names(exons.rbm17)
## freqs <- c(100, 200, 300)
## names(freqs) <- nms
## p.splice1 <- qplot(exons.rbm17)
## ## when turning on frequency 
## ## p.splice <- qplot(exons.rbm17, freq = freqs, show.label = TRUE, label.type = "count",
## ##       scale.size = c(1, 5), label.size = 3)
## p.splice2 <- qplot(exons.rbm17, freq = freqs, show.label = TRUE, offset = 0.05,
##                    label.type = "count")


###################################################
### code chunk number 14: iranges-plot (eval = FALSE)
###################################################
## ir <- ranges(gr[seqnames(gr) == "chr1"])[1:40]
## p1 <- qplot(ir) + opts(title = "full")
## p2 <- qplot(ir, geom = "segment")+ opts(title = "segment")
## p3 <- qplot(ir, geom = "coverage.line")+ opts(title = "coverage.line")
## p4 <- qplot(ir, geom = "coverage.polygon")+ opts(title = "coverage.polygon")
## grid.arrange(p1, p2, p3, p4, ncol = 2)


###################################################
### code chunk number 15: Rlelist-qplot
###################################################
set.seed(1)
lambda <- c(rep(0.001, 4500), seq(0.001, 10, length = 500), 
            seq(10, 0.001, length = 500))
xVector <- rpois(1e4, lambda)
xRle <- Rle(xVector)
xRleList <- RleList(xRle, 2L * xRle)
p1 <- qplot(xRle) + opts(title = "raw(point)")
p2 <- qplot(xRle, geom = "line")+ opts(title = "raw(line)")
p3 <- qplot(xRle, geom = "segment")+ opts(title = "raw(segment)")
p4 <- qplot(xRle, type = "viewMaxs", geom = "line", lower = 5)+
   opts(title = "viewMax(line)")
p5 <- qplot(xRle, type = "viewMins", geom = "line",lower = 5)+
   opts(title = "viewMins(line)")
p6 <- qplot(xRle, type = "viewMeans", geom = "line",lower = 5) +
   opts(title = "viewMeans(line)")
p7 <- qplot(xRle, type = "viewSums",geom = "line", lower = 5)+
   opts(title = "viewSums(line)")

## more control like this
## qplot(xRle, type = "viewSums", lower = 5, size = I(10), color = I("red"),
##      alpha = y)


###################################################
### code chunk number 16: rletrakcs
###################################################
grid.arrange(p1, p2, p3, p4, p5, p6, p7, ncol = 2)


###################################################
### code chunk number 17: rlelist-multiplesample
###################################################
p <- qplot(xRleList, geom = "line", facetByRow = TRUE)
print(p)


###################################################
### code chunk number 18: gapped-plot
###################################################
library(Rsamtools)
bamfile <- system.file("extdata", "SRR027894subRBM17.bam", package="biovizBase")
## need to set use.names = TRUE
ga <- readBamGappedAlignments(bamfile, use.names = TRUE) 
p1 <- qplot(ga)
p2 <- qplot(ga, show.junction = TRUE)
p3 <- qplot(ga, geom = "full")
grid.arrange(p1, p2, p3, ncol = 1)


###################################################
### code chunk number 19: bam-plot (eval = FALSE)
###################################################
## p1 <- qplot(bamfile, which = genesymbol["RBM17"], geom = "gapped.pair")
## p2 <- qplot(bamfile, which = genesymbol["RBM17"], geom = "gapped.pair",
##             show.junction = TRUE) 
## p3 <- qplot(bamfile, which = genesymbol["RBM17"], geom = "full")
## p4 <- qplot(bamfile, which = genesymbol["RBM17"], geom = "coverage.line") 
## p5 <- qplot(bamfile, which = genesymbol["RBM17"], geom = "coverage.polygon")
## ## p6 <- qplot(bamfile, which = genesymbol["RBM17"], bsgenome = Hsapiens,
## ##       geom = "mismatch.summary") 
## tracks(p1, p2, p3, p4, p5)


###################################################
### code chunk number 20: txdb-plot (eval = FALSE)
###################################################
## library(ggbio)
## data(genesymbol)
## pdf("intro-txdb-plot")
## p.full <- qplot(txdb, geom = "full", which = genesymbol["RBM17"])
## p.single <- qplot(txdb, geom = "single", which = genesymbol["RBM17"])
## p.tx <- qplot(txdb, geom = "tx", which = genesymbol["RBM17"])
## tracks(p.full, p.single, p.tx, heights = c(400, 100, 400))
## dev.off()


###################################################
### code chunk number 21: BSgenome-tracks
###################################################
library(BSgenome.Hsapiens.UCSC.hg19)
gr <- GRanges("chr1", IRanges(5e7, 5e7+50))
p1 <- qplot(Hsapiens, name = gr, geom = "text") +
  theme_bw()
p2 <- qplot(Hsapiens, name = gr, geom = "point") +
    theme_bw()
p3 <- qplot(Hsapiens, name = gr, geom = "segment") +
    theme_bw()
p4 <- qplot(Hsapiens, name = gr, geom = "rectangle") +
    theme_bw()
tracks(p1, p2, p3, p4)


###################################################
### code chunk number 22: plotOverview-cyto
###################################################
data(hg19IdeogramCyto)
## make shorter and clean labels
old.chrs <- seqnames(seqinfo(hg19IdeogramCyto))
new.chrs <- gsub("chr", "", old.chrs)
## lst <- as.list(new.chrs)
names(new.chrs) <- old.chrs
new.ideo <- renameSeqlevels(hg19IdeogramCyto, new.chrs)
p <- plotOverview(new.ideo, cytoband = TRUE)
print(p)


###################################################
### code chunk number 23: plotOverview-nocyto
###################################################
p <- plotOverview(new.ideo, cytoband = FALSE)
print(p)


###################################################
### code chunk number 24: plotOverview-nocyto-darned
###################################################
data(darned_hg19_subset500)
## rename 
old.chrs <- seqnames(seqinfo(darned_hg19_subset500))
new.chrs <- gsub("chr", "", old.chrs)
names(new.chrs) <- old.chrs
new.darned <- renameSeqlevels(darned_hg19_subset500, new.chrs)
p <- p + geom_hotregion(new.darned)
print(p)


###################################################
### code chunk number 25: plotOverview-nocyto-darned-exon
###################################################
p <- plotOverview(new.ideo, cytoband = FALSE)
p <- p + geom_hotregion(new.darned, aes(color = exReg))
print(p)


###################################################
### code chunk number 26: plotSingleChrom (eval = FALSE)
###################################################
## p <- plotSingleChrom(hg19IdeogramCyto, subchr = "chr1")
## p.zoom <- plotSingleChrom(hg19IdeogramCyto, subchr = "chr1",
##                 zoom.region = c(1e8, 1.5e8))
## 
## ## you could also run the following code to make nice adjusted 
## ## overview
## ## vp1 <- viewport(width = 1, height = 0.14)
## ## p <- plotSingleChrom(hg19IdeogramCyto, subchr = "chr1")
## ## print(p, vp = vp1)


###################################################
### code chunk number 27: manhattan-plot
###################################################
data(hg19Ideogram)
chrs <- as.character(levels(seqnames(hg19IdeogramCyto)))
seqlths <- seqlengths(hg19Ideogram)[chrs]
set.seed(1)
nchr <- length(chrs)
nsnps <- 500
gr.snp <- GRanges(rep(chrs,each=nsnps),
             IRanges(start = do.call(c, lapply(chrs, function(chr){
               N <- seqlths[chr]
               runif(nsnps,1,N)
             })), width = 1),
             SNP=sapply(1:(nchr*nsnps), function(x) paste("rs",x,sep='')),
             pvalue =  -log10(runif(nchr*nsnps)),
                  group = sample(c("Normal", "Tumor"), size = nchr*nsnps,
                    replace = TRUE)
             )

## processing the name, to make it shorter
nms <- seqnames(seqinfo(gr.snp))
nms.new <- gsub("chr", "", nms)
names(nms.new) <- nms
gr.snp <- renameSeqlevels(gr.snp, nms.new)
p1 <- plotGrandLinear(gr.snp, y = pvalue, geom = "point")
p2 <- plotGrandLinear(gr.snp, y = pvalue, geom = "point", color.type = "identity",
                color = I("black"), size = I(1))
## remove ticks and laebels and set cutoff
p3 <- plotGrandLinear(gr.snp, y = pvalue, geom = "point", cutoff = 4) +
            opts(axis.text.x = theme_blank(),
                 axis.ticks=theme_blank())
grid.arrange(p1, p2, p3, ncol = 1)


###################################################
### code chunk number 28: overview-facet
###################################################
p <- plotGrandLinear(gr.snp, y = pvalue, geom = "point",
                facet = group ~ .,  color.type = "twocolor")
print(p)


###################################################
### code chunk number 29: plotSpliceSum (eval = FALSE)
###################################################
## p1 <- plotSpliceSum(bamfile, exons.rbm17)
## ## all more control
## p2 <- plotSpliceSum(bamfile, exons.rbm17, weighted = FALSE, offset = 0.01)
## p3 <- plotSpliceSum(bamfile, txdb, which = genesymbol["RBM17"],
##                    show.label = TRUE,
##                    label.type = "count")
## print(p1)


###################################################
### code chunk number 30: frag-length (eval = FALSE)
###################################################
## model <- exonsBy(txdb, by = "tx")
## model.new <- subsetByOverlaps(model, genesymbol["RBM17"])
## exons.rbm17 <- subsetByOverlaps(exons(txdb), genesymbol["RBM17"])
## exons.new <- reduce(exons.rbm17)
## ## the graphic shows here is using different bamfile from the one
## ## in ext, it contains more pair-end reads.
## plotFragLength(bamfile, exons.new, geom = c("point","segment"))
## plotFragLength(bamfile, exons.new, geom = c("point","segment"), 
##                      annotation = FALSE)
## 
## plotFragLength(bamfile, exons.new, geom = c("point","segment"), 
##                      type = "cut", gap.ratio = 0.001)


###################################################
### code chunk number 31: mismatch-plot
###################################################
gr <- GRanges("chr10", IRanges(6134000, 6135000))
test <- pileupAsGRanges(bamfile, region = gr)
test.match <- pileupGRangesAsVariantTable(test, Hsapiens)
## use plotMismatchSum directly
pmis1 <- plotMismatchSum(test.match, show.coverage = FALSE)
pmis2 <- plotMismatchSum(test.match, show.coverage = TRUE)
grid.arrange(pmis1, pmis2, ncol = 1)


###################################################
### code chunk number 32: mismatch2 (eval = FALSE)
###################################################
## p <- qplot(test.match, geom = "mismatch.summary")
## library(Rsamtools)
## ## for character
## p <- qplot(bamfile, which = gr, bsgenome = Hsapiens,
##       geom = "mismatch.summary", show.coverage = TRUE)
## ## for BamFile
## p <- qplot(BamFile(bamfile), which = gr, bsgenome = Hsapiens,
##       geom = "mismatch.summary", show.coverage = TRUE)


###################################################
### code chunk number 33: plotRangesLinkedToData (eval = FALSE)
###################################################
## exons <- exons(txdb)
## exon17 <- subsetByOverlaps(exons, genesymbol["RBM17"])
## ## reduce to make sure there is no overlap
## ## just for example
## exon.new <- reduce(exon17)
## p <- qplot(model.new)
## ## simulated data
## values(exon.new)$sample1 <- rnorm(length(exon.new), 10, 3)
## values(exon.new)$sample2 <- rnorm(length(exon.new), 10, 10)
## values(exon.new)$score <- rnorm(length(exon.new))
## plotRangesLinkedToData(exon.new, stat.col = c("sample1", "sample2"))
## ## show as figure, no annotation
## ## the same as
## ## plotRangesLinkedToData(exon.new, stat.col = 1:2)
## ## show as figure, adding annotation


###################################################
### code chunk number 34: link2 (eval = FALSE)
###################################################
## plotRangesLinkedToData(exon.new, stat.col = 1:2, annotation = list(p))


###################################################
### code chunk number 35: sessionInfo
###################################################
sessionInfo()