## ----style-knitr, eval=TRUE, echo=FALSE, results="asis"------------------
BiocStyle::latex()

## ----setup, include=FALSE, cache=FALSE, eval = TRUE, echo = FALSE--------
library(knitr)
opts_chunk$set(fig.path='./figures/ggbio-',
               fig.align='center', fig.show='asis',
               eval = TRUE,
               fig.width = 4.5,
               fig.height = 4.5,
               tidy = FALSE,
               message = FALSE,
               warning = FALSE)
options(replace.assign=TRUE,width=80)

## ----citation------------------------------------------------------------
citation("ggbio")

## ------------------------------------------------------------------------

## ------------------------------------------------------------------------

## ------------------------------------------------------------------------

## ----fig.width = 5.5, fig.height = 1.5-----------------------------------
library(ggbio)
p.ideo <- Ideogram(genome = "hg19")
p.ideo
library(GenomicRanges)
## special highlights instead of zoomin!
p.ideo + xlim(GRanges("chr2", IRanges(1e8, 1e8+10000000)))

## ------------------------------------------------------------------------
library(ggbio)
library(Homo.sapiens)
class(Homo.sapiens)
##
data(genesymbol, package = "biovizBase")
wh <- genesymbol[c("BRCA1", "NBR1")]
wh <- range(wh, ignore.strand = TRUE)


p.txdb <- autoplot(Homo.sapiens, which  = wh)
p.txdb

autoplot(Homo.sapiens, which  = wh, label.color = "black", color = "brown",
          fill = "brown")

## ------------------------------------------------------------------------
autoplot(Homo.sapiens, which  = wh, gap.geom = "chevron")

## ----fig.width = 5.5, fig.height = 1.5-----------------------------------
autoplot(Homo.sapiens, which = wh, stat = "reduce")

## ------------------------------------------------------------------------
columns(Homo.sapiens)
autoplot(Homo.sapiens, which  = wh, columns = c("TXNAME", "GO"), names.expr = "TXNAME::GO")

## ------------------------------------------------------------------------
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
autoplot(txdb, which = wh)

## ------------------------------------------------------------------------
library(biovizBase)
gr.txdb <- crunch(txdb, which = wh)
## change column to 'model'
colnames(values(gr.txdb))[4] <- "model"
grl <- split(gr.txdb, gr.txdb$tx_id)
## fake some randome names
names(grl) <- sample(LETTERS, size = length(grl), replace = TRUE)
grl

## ------------------------------------------------------------------------
autoplot(grl, aes(type = model))
ggplot() + geom_alignment(grl, type = "model")

## ----fig.width = 5.5, fig.height = 1.5-----------------------------------
library(BSgenome.Hsapiens.UCSC.hg19)
bg <- BSgenome.Hsapiens.UCSC.hg19
p.bg <- autoplot(bg, which = wh)
## no geom
p.bg
## segment
p.bg + zoom(1/100)
## rectangle
p.bg + zoom(1/1000)
## text
p.bg + zoom(1/2500)

## ----eval = FALSE--------------------------------------------------------
## library(BSgenome.Hsapiens.UCSC.hg19)
## bg <- BSgenome.Hsapiens.UCSC.hg19
## ## force to use geom 'segment' at this level
## autoplot(bg, which = resize(wh, width = width(wh)/2000), geom = "segment")

## ------------------------------------------------------------------------
fl.bam <- system.file("extdata", "wg-brca1.sorted.bam", package  = "biovizBase")
wh <- keepSeqlevels(wh, "chr17")
autoplot(fl.bam, which = wh)

## ------------------------------------------------------------------------
fl.bam <- system.file("extdata", "wg-brca1.sorted.bam", package  = "biovizBase")
wh <- keepSeqlevels(wh, "chr17")
autoplot(fl.bam, which = resize(wh, width = width(wh)/10), geom = "gapped.pair")

## ------------------------------------------------------------------------
library(BSgenome.Hsapiens.UCSC.hg19)
bg <- BSgenome.Hsapiens.UCSC.hg19
p.mis <- autoplot(fl.bam, bsgenome = bg, which = wh, stat = "mismatch")
p.mis

## ------------------------------------------------------------------------
autoplot(fl.bam, method = "estimate")
autoplot(fl.bam, method = "estimate", which = paste0("chr", 17:18), aes(y = log(..coverage..)))

## ------------------------------------------------------------------------
library(VariantAnnotation)
fl.vcf <- system.file("extdata", "17-1409-CEU-brca1.vcf.bgz", package="biovizBase")
vcf <- readVcf(fl.vcf, "hg19")
vr <- as(vcf[, 1:3], "VRanges")
vr <- renameSeqlevels(vr, value = c("17" = "chr17"))
## small region contains data
gr17 <- GRanges("chr17", IRanges(41234400, 41234530))
p.vr <- autoplot(vr, which = wh)
## none geom
p.vr
## rect geom
p.vr + xlim(gr17)
## text geom
p.vr + xlim(gr17) + zoom()

## ----eval = FALSE--------------------------------------------------------
## autoplot(vr, which = wh, geom = "rect", arrow = FALSE)

## ----fig.width = 8, fig.height = 5.5-------------------------------------
## tks <- tracks(p.ideo, mismatch = p.mis, dbSNP = p.vr, ref = p.bs, gene = p.txdb)
## tks <- tracks(fl.bam, fl.vcf, bs, Homo.sapiens) ## default ideo = FALSE, turned on
## tks <- tracks(fl.bam, fl.vcf, bs, Homo.sapiens, ideo = TRUE)
## tks + xlim(gr17)
gr17 <- GRanges("chr17", IRanges(41234415, 41234569))
tks <- tracks(p.ideo, mismatch = p.mis, dbSNP = p.vr, ref = p.bg, gene = p.txdb,
              heights = c(2, 3, 3, 1, 4)) + xlim(gr17) + theme_tracks_sunset()
tks

## ----fig.width = 8, fig.height = 5.5-------------------------------------
## zoom in
tks + zoom()

## ----eval = FALSE--------------------------------------------------------
## ## zoom in with scale
## p.txdb + zoom(1/8)
## ## zoom out
## p.txdb + zoom(2)
## ## next view page
## p.txdb + nextView()
## ## previous view page
## p.txdb + prevView()

## ----processing----------------------------------------------------------
data("CRC", package  = "biovizBase")

## ------------------------------------------------------------------------
head(hg19sub)
autoplot(hg19sub, layout = "circle", fill  = "gray70")

## ------------------------------------------------------------------------
p <- ggbio() + circle(hg19sub, geom = "ideo", fill = "gray70") +
    circle(hg19sub, geom = "scale", size = 2) +
  circle(hg19sub, geom = "text", aes(label = seqnames), vjust = 0, size = 3)
p

## ------------------------------------------------------------------------
p <- ggbio(trackWidth = 10, buffer = 0, radius = 10) + circle(hg19sub, geom = "ideo", fill = "gray70") +
    circle(hg19sub, geom = "scale", size = 2) +
  circle(hg19sub, geom = "text", aes(label = seqnames), vjust = 0, size = 3)
p

## ----lower-mut-track-----------------------------------------------------
head(mut.gr)
p <- ggbio() + circle(mut.gr, geom = "rect", color = "steelblue") +
    circle(hg19sub, geom = "ideo", fill = "gray70") +
    circle(hg19sub, geom = "scale", size = 2) +
  circle(hg19sub, geom = "text", aes(label = seqnames), vjust = 0, size = 3)
p

## ------------------------------------------------------------------------
head(crc.gr)

## ----subset-crc-1--------------------------------------------------------
gr.crc1 <- crc.gr[values(crc.gr)$individual == "CRC-1"]

## ----lower-point-track---------------------------------------------------
## manually specify radius
p <- p + circle(gr.crc1, geom = "point", aes(y = score, size = tumreads),
                color = "red", grid = TRUE, radius = 30) + scale_size(range = c(1, 2.5))
p

## ----lower-link-track----------------------------------------------------
## specify radius manually
p <- p + circle(gr.crc1, geom = "link", linked.to = "to.gr", aes(color = rearrangements),
                       radius = 23)
p

## ------------------------------------------------------------------------
p <- ggbio() +
   circle(gr.crc1, geom = "link", linked.to = "to.gr", aes(color = rearrangements)) +
  circle(gr.crc1, geom = "point", aes(y = score, size = tumreads),
                color = "red", grid = TRUE) + scale_size(range = c(1, 2.5)) +
  circle(mut.gr, geom = "rect", color = "steelblue") +
    circle(hg19sub, geom = "ideo", fill = "gray70") +
    circle(hg19sub, geom = "scale", size = 2) +
  circle(hg19sub, geom = "text", aes(label = seqnames), vjust = 0, size = 3)
p

## ----arrangement---------------------------------------------------------
grl <- split(crc.gr, values(crc.gr)$individual)
## need "unit", load grid
library(grid)
crc.lst <- lapply(grl, function(gr.cur){
  print(unique(as.character(values(gr.cur)$individual)))
  cols <- RColorBrewer::brewer.pal(3, "Set2")[2:1]
  names(cols) <- c("interchromosomal", "intrachromosomal")
  p <- ggbio() + circle(gr.cur, geom = "link", linked.to = "to.gr",
                         aes(color = rearrangements)) +
                  circle(hg19sub, geom = "ideo",
                         color = "gray70", fill = "gray70") +
                  scale_color_manual(values = cols)  +
                  labs(title = (unique(values(gr.cur)$individual))) +
                  theme(plot.margin = unit(rep(0, 4), "lines"))
})

## ----simple-wrapper, fig.width = 7, fig.height = 5-----------------------

arrangeGrobByParsingLegend(crc.lst, widths = c(4, 1), legend.idx = 1, ncol = 3)

## ----simul_snp-----------------------------------------------------------
snp <- read.table(system.file("extdata", "plink.assoc.sub.txt", package = "biovizBase"),
                  header = TRUE)
require(biovizBase)
gr.snp <- transformDfToGr(snp, seqnames = "CHR", start = "BP", width = 1)
head(gr.snp)
## change the seqname order
require(GenomicRanges)
gr.snp <- keepSeqlevels(gr.snp, as.character(1:22))
seqlengths(gr.snp)
## need to assign seqlengths
data(ideoCyto, package = "biovizBase")
seqlengths(gr.snp) <- as.numeric(seqlengths(ideoCyto$hg18)[1:22])
## remove missing
gr.snp <- gr.snp[!is.na(gr.snp$P)]
## transform pvalue
values(gr.snp)$pvalue <- -log10(values(gr.snp)$P)
head(gr.snp)
## done

## ----line, fig.height = 4, fig.width = 6---------------------------------
autoplot(gr.snp, geom = "point", coord = "genome", aes(y = pvalue))

## ----morecolor2, fig.height = 4, fig.width = 6---------------------------
plotGrandLinear(gr.snp, aes(y = pvalue), color = c("#7fc97f", "#fdc086"))

## ----fig.height = 4, fig.width = 6---------------------------------------
plotGrandLinear(gr.snp, aes(y = pvalue), color = c("#7fc97f", "#fdc086"),
                cutoff = 3, cutoff.color = "blue", cutoff.size = 0.2)


## ----fig.height = 4, fig.width = 6---------------------------------------

plotGrandLinear(gr.snp, aes(y = pvalue, color = OR), spaceline = TRUE, legend = TRUE)

## ----fig.height = 4, fig.width = 6---------------------------------------
gro <- GRanges(c("1", "11"), IRanges(c(100, 2e6), width = 5e7))
names(gro) <- c("group1", "group2")
plotGrandLinear(gr.snp, aes(y = pvalue), highlight.gr = gro)

## ------------------------------------------------------------------------
data(ideoCyto, package = "biovizBase")
autoplot(seqinfo(ideoCyto$hg19), layout = "karyogram")

## ----fig.width = 5, fig.height = 4---------------------------------------
## turn on cytoband if it exists
biovizBase::isIdeogram(ideoCyto$hg19)
autoplot(ideoCyto$hg19, layout = "karyogram", cytoband = TRUE)

## ----load-RNAediting-----------------------------------------------------
data(darned_hg19_subset500, package = "biovizBase")
dn <- darned_hg19_subset500
library(GenomicRanges)
seqlengths(dn)
## add seqlengths
## we have seqlegnths information in another data set
seqlengths(dn) <- seqlengths(ideoCyto$hg19)[names(seqlengths(dn))]
## then we change order
dn <- keepSeqlevels(dn, paste0("chr", c(1:22, "X")))
seqlengths(dn)
autoplot(dn, layout = "karyogram")

## ----load-RNAediting-color-----------------------------------------------
## since default is geom rectangle, even though it's looks like segment
## we still use both fill/color to map colors
autoplot(dn, layout = "karyogram", aes(color = exReg, fill = exReg))

## ------------------------------------------------------------------------
## since default is geom rectangle, even though it's looks like segment
## we still use both fill/color to map colors
autoplot(dn, layout = "karyogram", aes(color = exReg, fill = exReg), alpha  = 0.5) +
  scale_color_discrete(na.value = "brown")

## ----fig.width = 5, fig.height = 4---------------------------------------
## let's remove the NA value
dn.nona <- dn[!is.na(dn$exReg)]
## compute levels based on categories
dn.nona$levels <- as.numeric(factor(dn.nona$exReg))
## do a trcik show them at different height
p.ylim <- autoplot(dn.nona, layout = "karyogram", aes(color = exReg, fill = exReg,
                                       ymin = (levels - 1) * 10/3,
                                       ymax = levels * 10 /3))

## ------------------------------------------------------------------------
## prepare the data
dn3 <- dn.nona[dn.nona$exReg == '3']
dn5 <- dn.nona[dn.nona$exReg == '5']
dnC <- dn.nona[dn.nona$exReg == 'C']
dn.na <- dn[is.na(dn$exReg)]
## now we have 4 different data sets
autoplot(seqinfo(dn3), layout = "karyogram") +
  layout_karyogram(data = dn3, geom = "rect", ylim = c(0, 10/3), color = "#7fc97f") +
  layout_karyogram(data = dn5, geom = "rect", ylim = c(10/3, 10/3*2), color = "#beaed4") +
  layout_karyogram(data = dnC, geom = "rect", ylim = c(10/3*2, 10), color = "#fdc086") +
  layout_karyogram(data = dn.na, geom = "rect", ylim = c(10, 10/3*4), color = "brown")

## ----edit-space, fig.width = 5, fig.height = 4---------------------------
dn$pvalue <- runif(length(dn)) * 10
p <- autoplot(seqinfo(dn)) + layout_karyogram(dn, aes(x = start, y = pvalue),
                     geom = "point", color = "#fdc086")
p

## ----fig.width = 5, fig.height = 4---------------------------------------
p.ylim + facet_wrap(~seqnames)

## ----fig.width = 8-------------------------------------------------------
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(ggbio)
data(genesymbol, package = "biovizBase")
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
model <- exonsBy(txdb, by = "tx")
model17 <- subsetByOverlaps(model, genesymbol["RBM17"])
exons <- exons(txdb)
exon17 <- subsetByOverlaps(exons, genesymbol["RBM17"])
## reduce to make sure there is no overlap
## just for example
exon.new <- reduce(exon17)
## suppose
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))
values(exon.new)$significant <- sample(c(TRUE,FALSE), size = length(exon.new),replace = TRUE)
## data ready
exon.new

## ------------------------------------------------------------------------
p17 <- autoplot(txdb, genesymbol["RBM17"])
plotRangesLinkedToData(exon.new, stat.y = c("sample1", "sample2"), annotation = list(p17))

## ------------------------------------------------------------------------
p.txdb
p.txdb + theme_alignment()
p.txdb + theme_clear()
p.txdb + theme_null()

## ------------------------------------------------------------------------
library(GenomicRanges)
set.seed(1)
N <- 100
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),
              sample = sample(c("Normal", "Tumor"),
                size = N, replace = TRUE),
              pair = sample(letters, size = N,
                replace = TRUE))
seqlengths(gr) <- c(400, 1000, 500)
autoplot(gr)
autoplot(gr) + theme_genome()

## ------------------------------------------------------------------------
theme_tracks_sunset

## ----session-info--------------------------------------------------------
sessionInfo()