##----------------------------------------------------------------------------##
## For class "CoverageView"
##----------------------------------------------------------------------------##
## coverage/alignment, separate model?
## make properties
gparslst <- list(view = "character",
lower = "numeric",
binNum = "PositiveInteger",
zoomLevel.cur = "numeric",
zoomLevel = "numeric",
geom = .GeomName("CoverageView"))
setPars("CoverageView", gparslst, methods = list(
widget = function(visible = "binNum"){
ControlPanel(.self, visible = visible)
}
))
CoverageView.gen <- setRefClass("CoverageView",
contains = c("QtVisnabView", "LinearView"),
fields=list(
track = "list",
file = "character",
BSgenome = "BSgenomeORNULL",
gr.v = "GRanges",
res = "data.frame",
cached = "logical",
cachedYlim = "numeric",
cachedSeqnames = "character",
cachedCoverage = "GRanges"
))
##----------------------------------------------------------------------##
## "CoverageView" constructor
##----------------------------------------------------------------------##
CoverageView <- function(file,
seqname,
BSgenome = NULL,
geom = c("total","mismatch","pairend","elength"),
rescale = c("none", "geometry","transform"),
viewname = "Coverage",
lower = 10L,
binNum = 50L,
hint = FALSE, #not implemented
...){
## get params
geom <- match.arg(geom)
geom <- new("CoverageViewGeomSingleEnum", geom)
rescale <- match.arg(rescale)
rescale <- new("RescaleSingleEnum", rescale)
tooltips <- "not implemented"
hd <- scanBamHeader(file)
.levels <- unique(names(hd[[1]]$targets))
if(missing(seqname))
seqname <- sort(names(hd[[1]]$targets))[1]
seqlength <- hd[[1]]$targets
xlimZoom <- c(1, seqlength[seqname])
viewrange <- MutableGRanges(factor(seqname, levels = sort(.levels)),
IRanges(1, seqlength[seqname]),
seqlengths = seqlength)
## seqlengths(viewrange) <- seqlength
## viewrange$seqinfo@seqlengths <- as.integer(seqlength)
pars <- GraphicPars(geom = geom,
lower = lower, binNum = PositiveInteger(binNum),
view = "CoverageView")
md <- IModeGroupWidgetQt(scaleMode = ScaleMode(zoomMode = "Horizontal"))
message("Estimating coverage...")
covg <- estimateCoverage(BamFile(file))
obj <- CoverageView.gen$new(file = file,
xlimZoom = xlimZoom,
BSgenome = BSgenome, mode = md,
rescale = rescale,
eventTrace = new("EventTrace"),
viewrange = viewrange,
pars = pars,
cached = FALSE,
cachedCoverage = covg
)
obj$createView()
obj$regSignal()
obj
}
##---------------------------------------------------------##
## print method
##---------------------------------------------------------##
CoverageView.gen$methods(createView = function(){
## TODO: for all??
cachedSeqnames <<- seqname <- as.character(seqnames(viewrange))
setDislayWidgets()
setBgColor()
hd <- scanBamHeader(file)
xlim <<- c(1, seqlengths(viewrange)[seqname])
pars$zoomLevel <<- c(1e6, 5e4, 5e4-1, 1e3, 0)
pars$zoomLevel.cur <<- 1
cur.gr <- cachedCoverage[seqnames(cachedCoverage) == seqname]
xpos <- (start(cur.gr) + end(cur.gr))/2
ypos <- cur.gr$score
idx <- order(xpos)
xpos <- xpos[idx]
ypos <- ypos[idx]
res <<- data.frame(xpos = xpos, ypos = ypos)
ylim <<- expand_range(c(0, max(ypos)), mul = 0.05)
gr.v <<- GRanges(seqnames(viewrange)@values, IRanges(start = 1, width = 1) )
pfunCov <- function(layer,painter,exposed){
xlimZoomChanged$block()
xlimZoom <<- as.matrix(exposed)[,1]
## viewrange$ranges <<- xlimZoom
if(!eventTrace$selfSignal){
## viewrange$changed$connect(function)
signal(viewrange)$unblock()
ranges(viewrange) <<- IRanges(xlimZoom[1] , xlimZoom[2])
}
if(eventTrace$selfSignal){
signal(viewrange)$block()
ranges(viewrange) <<- IRanges(xlimZoom[1] , xlimZoom[2])
}
xlimZoomChanged$unblock()
xlimZoom <<- as.matrix(exposed)[,1]
##======================================================================
## geom == "total"
##======================================================================
if(pars$geom == "total"){
## level 1
if(diff(xlimZoom)>pars$zoomLevel[1]){
if(pars$zoomLevel.cur != 1){
cur.gr <- cachedCoverage[seqnames(cachedCoverage) == seqname]
xpos <- (start(cur.gr) + end(cur.gr))/2
ypos <- cur.gr$score
idx <- order(xpos)
xpos <- xpos[idx]
ypos <- ypos[idx]
res <<- data.frame(xpos = xpos, ypos = ypos)
ylim <<- expand_range(c(0, max(ypos)), mul = 0.05)
}
.self$paintCovSeg(painter, res)
pars$zoomLevel.cur <<- 1
}
## level 2
if(diff(xlimZoom)<=pars$zoomLevel[1] &
diff(xlimZoom)>pars$zoomLevel[2]){
if(pars$zoomLevel.cur != 2){
cur.gr <- cachedCoverage[seqnames(cachedCoverage) == seqname]
xpos <- (start(cur.gr) + end(cur.gr))/2
ypos <- cur.gr$score
idx <- order(xpos)
xpos <- xpos[idx]
ypos <- ypos[idx]
res <<- data.frame(xpos = xpos, ypos = ypos)
ylim <<- expand_range(c(0, max(ypos)), mul = 0.05)
}
.self$paintCovHist(painter, res, bin = 1e4)
pars$zoomLevel.cur <<- 2
}
## level 3
## if(diff(xlimZoom)<=pars$zoomLevel[2] &
## diff(xlimZoom)>pars$zoomLevel[3]){
## if(pars$zoomLevel.cur != 3){
## .diff <- diff(xlimZoom)
## message("computing bam count...")
## ## nbin <- (end(range(gr.v)) - start(range(gr.v)))/
## ## (diff(xlimZoom)/pars$binNum)
## nbin <- 100
## sts <- seq(xlimZoom[1]-0.25 * .diff, xlimZoom[2] + 0.25 * .diff,
## length.out = nbin)
## gr.v <<- GRanges(seqnames(viewrange)@values,
## IRanges(start = sts, width = diff(sts)[1]))
## res <<- countBam(file, param = ScanBamParam(which = gr.v))
## }else{
## .diff <- diff(xlimZoom)
## if(xlimZoom[2] >= max(end(gr.v)) | xlimZoom[1] <= min(start(gr.v))){
## message("clear cache... recomputing bam count...")
## nbin <- 100
## sts <- seq(xlimZoom[1]- 0.25*.diff, xlimZoom[2] + 0.25*.diff,
## length.out = nbin)
## gr.v <<- GRanges(seqnames(viewrange)@values,
## IRanges(start = sts, width = diff(sts)[1]))
## res <<- countBam(file, param = ScanBamParam(which = gr.v))
## }
## }
## .self$paintCovRect(painter, res)
## pars$zoomLevel.cur <<- 3
## }
## level 4
if(diff(xlimZoom)<=pars$zoomLevel[3] &
diff(xlimZoom)>pars$zoomLevel[4]){
if(pars$zoomLevel.cur != 4){
.diff <- diff(xlimZoom)
ir <- ranges(viewrange)
gr.v <<- GRanges(seqnames(viewrange),
IRanges(start(ir) - .diff, end(ir) + .diff))
suppressMessages(temp <- biovizBase:::fetch(BamFile(file), which = gr.v, type = "raw"))
temp <- keepSeqlevels(temp, as.character(unique(seqnames(temp))))
seqlengths(temp) <- end(gr.v)
covs <- coverage(temp)[[1]]
ypos <- as.numeric(covs[start(gr.v):end(gr.v)])
xpos <- start(gr.v):end(gr.v)
res <<- data.frame(xpos = xpos, ypos = ypos)
}else{
.diff <- diff(xlimZoom)
if(xlimZoom[2] >= max(end(gr.v)) | xlimZoom[1] <= min(start(gr.v))){
ir <- ranges(viewrange)
gr.v <<- GRanges(seqnames(viewrange),
IRanges(start(ir) - .diff, end(ir) + .diff))
suppressMessages(temp <- biovizBase:::fetch(BamFile(file), which = gr.v, type = "raw"))
temp <- keepSeqlevels(temp, as.character(unique(seqnames(temp))))
seqlengths(temp) <- end(gr.v)
covs <- coverage(temp)[[1]]
ypos <- as.numeric(covs[start(gr.v):end(gr.v)])
xpos <- start(gr.v):end(gr.v)
res <<- data.frame(xpos = xpos, ypos = ypos)
}
}
.self$paintCovHist(painter, res)
pars$zoomLevel.cur <<- 4
}
## level 5
if(diff(xlimZoom)<=pars$zoomLevel[4] &
diff(xlimZoom)>pars$zoomLevel[5]){
if(pars$zoomLevel.cur != 5){
.diff <- diff(xlimZoom)
message("computing bam stepping...")
ir <- ranges(viewrange)
gr.v <<- GRanges(seqnames(viewrange),
IRanges(start(ir) - .diff, end(ir) + .diff))
suppressMessages(temp <-
biovizBase:::fetch(BamFile(file), which = gr.v, type = "raw"))
temp$stepping <- disjointBins(ranges(temp))
strs <- as.character(strand(temp))
temp$.color <- as.character(biovizBase::getBioColor("STRAND")[strs])
res <<- as.data.frame(temp)
message("done")
}else{
.diff <- diff(xlimZoom)
if(xlimZoom[2] >= max(end(gr.v)) | xlimZoom[1] <= min(start(gr.v))){
message("reload...")
ir <- ranges(viewrange)
gr.v <<- GRanges(seqnames(viewrange),
IRanges(start(ir) - .diff, end(ir) + .diff))
suppressMessages(temp <-
biovizBase:::fetch(BamFile(file), which = gr.v, type = "raw"))
strs <- as.character(strand(temp))
temp$.color <- as.character(biovizBase::getBioColor("STRAND")[strs])
temp$stepping <- disjointBins(ranges(temp))
res <<- as.data.frame(temp)
message("done")
}
}
.self$paintStep(painter, res)
pars$zoomLevel.cur <<- 5
}
}
##======================================================================
## geom == "mismatch"
##======================================================================
if(pars$geom == "mismatch"){
## level 1
if(diff(xlimZoom)>pars$zoomLevel[1]){
.self$paintCovSeg(painter, xpos, ypos)
}
if(diff(xlimZoom)<=pars$zoomLevel[1] &
diff(xlimZoom)>pars$zoomLevel[2]){
.self$paintCovRect(painter)
}
if(diff(xlimZoom)<=pars$zoomLevel[2] & diff(xlimZoom)>pars$zoomLevel[3]){
## this should require loading or associated BSgenome
## And support sequence logo(maybe seperate painter)
.self$paintCovMismatch(painter)
}}
if(pars$geom == "elength"){
if(diff(xlimZoom)> 1e5){
.self$paintCovSeg(painter, xpos, ypos)
}
if(diff(xlimZoom)<= 1e5 & diff(xlimZoom)>pars$zoomLevel[3]){
.self$paintFragLength(painter)
}
}
}
keyOutFun <- function(layer, event){
eventTrace$focusin <<- FALSE
}
hoverEnterFun <- function(layer, event){
eventTrace$focusin <<- TRUE
}
hoverLeaveFun <- function(layer, event){
## eventTrace$focusin <<- FALSE
eventTrace$focusin <<- TRUE
}
rootLayer[0,0] <<- qlayer(scene, paintFun=pfunCov,
row=row, col=col, rowSpan=rowSpan, colSpan=colSpan,
wheelFun = wheelEventZoom(),
keyPressFun = keyPressEventZoom(),
hoverEnterFun = hoverEnterFun,
focusOutFun = keyOutFun, hoverLeaveFun = hoverLeaveFun)
rootLayer[0,0]$setLimits(qrect(xlim, ylim))
ylimChanged$connect(function(){
rootLayer[0,0]$setLimits(qrect(xlim, ylim))
})
})
CoverageView.gen$methods(show = function(){
view$show()
})
setMethod("print","CoverageView",function(x,..){
x$show()
})
CoverageView.gen$methods(regSignal = function(){
pars$binNumChanged$connect(function(){
qupdate(scene)
})
xlimZoomChanged$connect(function(){
zoom_factor <- diff(xlimZoom)/seqlengths(viewrange)[cachedSeqnames]
## then scale view
view$resetTransform()
view$scale(1/zoom_factor, 1)
## then center viewr
pos.x <- mean(xlimZoom)
pos.y <- mean(ylim)
pos.scene <- as.numeric(rootLayer[0,0]$mapToScene(pos.x, pos.y))
view$centerOn(pos.scene[1], pos.scene[2])
ranges(viewrange) <<- IRanges(xlimZoom[1] , xlimZoom[2])
})
pars$geomChanged$connect(function(){
qupdate(scene)
})
viewrange <<- connect(viewrange, function(change){
if(rangeChanged(change)){
if(as.character(seqnames(viewrange)) != cachedSeqnames){
cachedSeqnames <<- as.character(seqnames(viewrange))
rootLayer[0,0]$close()
view$resetTransform()
createView()
regSignal()
}
}
})
theme$bgColorChanged$connect(function(){
bgcol <- pars$bgColor
bgalpha <- pars$alpha
qcol <- col2qcol(bgcol,bgalpha)
scene$setBackgroundBrush(qbrush(qcol))
})
})
## CoverageView.gen$methods(paintCovEst = function(painter, xpos, ypos){
## ## qdrawPolygon(painter,c(xpos[1], xpos, xpos[length(xpos)]),
## ## c(0, ypos, 0),stroke="gray40", fill = "gray40")
## qdrawSegment(painter,c(xpos[1], xpos, xpos[length(xpos)]),
## c(0, ypos, 0),stroke="gray40", fill = "gray40")
## })
CoverageView.gen$methods(paintCovSeg = function(painter, res){
qdrawSegment(painter,res$xpos,0, res$xpos, res$ypos,stroke="gray40")
ylim <<- expand_range(c(0, max(res$ypos)), mul = 0.05)
})
CoverageView.gen$methods(paintCovHist = function(painter, res, bin = 1){
qdrawRect(painter,res$xpos, 0, res$xpos + bin,
res$ypos,stroke="gray40", fill = "gray40")
ylim <<- expand_range(c(0, max(res$ypos)), mul = 0.05)
})
## CoverageView.gen$methods(paintCovPoly = function(painter, xpos, ypos){
## qdrawSegment(painter,xpos,0, xpos, ypos,stroke="gray40")
## })
## CoverageView.gen$methods(paintCovHist = function(painter, xpos, ypos){
## qdrawSegment(painter,xpos,0, xpos+, ypos,stroke="gray40")
## })
## CoverageView.gen$methods(paintCovCount = function(painter, xpos1, xpos2, ypos){
## qdrawRect(painter, xpos1, 0, xpos2, ypos, stroke = NA, fill = "gray40")
## })
CoverageView.gen$methods(paintCovRect = function(painter, res){
qdrawRect(painter, res$start, 0, res$end,
res$records, fill="gray40", stroke = "white")
ylim <<- expand_range(c(0, max(res$records)), mul = 0.05)
})
CoverageView.gen$methods(paintStep = function(painter, res){
## print(head(res$.color))
qdrawRect(painter, res$start, res$stepping -0.4, res$end,
res$stepping + 0.4, fill= as.character(res$.color),
stroke = as.character(res$.color))
ylim <<- expand_range(c(0, max(res$stepping)), mul = 0.05)
})
CoverageView.gen$methods(paintCovMismatch = function(painter){
if(is.null(BSgenome))
stop("Please specify the associated BSgenome object
if geom is set to mismatch")
gr <- GRanges(seqnames = seqnames(viewrange),
IRanges(xlimZoom[1], xlimZoom[2]))
lgr <- pileupAsGRanges(file, gr)
if(length(lgr)>0){
lgr <- pileupGRangesAsVariantTable(lgr, BSgenome) #too slow...
## !!! need to be removed
## idx <- values(lgr)$count/values(lgr)$depth>0.05 & !values(lgr)$match
idx <- values(lgr)$count/values(lgr)$depth>0.05 & !values(lgr)$match
idx <- sample(which(idx), size = as.integer(sum(idx)*0.9), replace = FALSE)
values(lgr)$count[idx] <- values(lgr)$count[idx]+as.integer(rnorm(1, mean = 10,
sd = 5))
values(lgr)$match[idx] <- TRUE
## group by match and mismatch
## idx <- values(lgr)$read != values(lgr)$ref
## lgr <- lgr[idx]
idx <- order(start(lgr), values(lgr)$match, values(lgr)$read)
## assumption: on the same chromosome
lgr <- lgr[idx]
eds <- unlist(lapply(split(values(lgr)$count, start(lgr)), function(x){
cumsum(x)
}))
sts <- unlist(lapply(split(values(lgr)$count, start(lgr)), function(x){
N <- length(x)
c(0,cumsum(x)[-N])
}))
ism <- unlist(lapply(split(values(lgr)$match, start(lgr)), function(x){
x
}))
rds <- unlist(lapply(split(values(lgr)$read, start(lgr)), function(x){
x
}))
## assign color
## cols <- genLegend(lgr, color = "read")$color
dnacol <- baseColor(IRanges:::safeExplode("ACTGN"))
miscols <- unlist(lapply(rds, function(rd){
dnacol[[rd]]
}))
miscols[ism] <- "gray40"
dnacol <- baseColor(IRanges:::safeExplode("ACTGN"))
cols <- unlist(lapply(rds, function(rd){
dnacol[[rd]]
}))
cols[ism] <- "gray40"
ylim <<- expand_range(range(log(eds + 1)), mul = 0.05)
## qdrawRect(painter, start(lgr), sts,
## start(lgr)+1, eds,
## stroke = NA, fill = miscols)
qdrawRect(painter, start(lgr), log(sts+1),
start(lgr)+1, log(eds+1),
stroke = NA, fill = cols)
}})
CoverageView.gen$methods(paintFragLength = function(painter){
gr <- GRanges(seqnames = seqnames(viewrange),
IRanges(xlimZoom[1], xlimZoom[2]))
pspan <- pspanGR(file, gr)$pspan
pspan <- subsetByOverlaps(ranges(pspan), ranges(gr), type = "within")
## draw point based gr
## THINK: MAKE IT GENERAL!
x <- (start(pspan)+end(pspan))/2
y <- width(pspan)
cir <- qglyphCircle(r = 1)
qdrawGlyph(painter, cir, x, y, stroke = NULL, fill = "black")
## qdrawCircle(painter, x, y, r = 1, stroke = NULL, fill = "black")
## fixed or not?
if(TRUE)
ylim <<- expand_range(c(min(y), max(y)), mul = 0.05)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.