#### Private Functions ####
###########################
#' .blankGr
#' @importFrom GenomicRanges makeGRangesFromDataFrame
.blankGr <- function(){
makeGRangesFromDataFrame(data.frame("chr"='chrZ', "start"=1, "end"=1))
}
#' .grepNA
#'
#' @param pattern pattern to grep
#' @param x vector to grep
.grepNA <- function(pattern, x){
idx <- grep(pattern, x)
if(length(idx) > 0) idx else NA
}
#' .getChrLength
#' @importFrom BSgenome.Hsapiens.UCSC.hg19 Hsapiens
#' @importFrom GenomicRanges makeGRangesFromDataFrame
#' @importFrom GenomeInfoDb seqlengths
#' @importFrom GenomicRanges end
.getChrLength <- function(){
chr.lengths = seqlengths(Hsapiens)[1:24]
chr.len.gr <- makeGRangesFromDataFrame(data.frame("chrom"=names(chr.lengths),
"loc.start"=rep(1, length(chr.lengths)),
"loc.end"=chr.lengths))
chr.len.gr$cum.end <- cumsum(as.numeric(end(chr.len.gr)))
chr.len.gr$cum.start <- chr.len.gr$cum.end - (end(chr.len.gr) -1)
chr.len.gr$cum.mid <- (chr.len.gr$cum.start + ((chr.len.gr$cum.end - chr.len.gr$cum.start)/2))
return(chr.len.gr)
}
#' .addCumPos
#'
#' @param dat Data containing $chrom, $pos, $loc.start, $loc.end
#' @param ref GRanges object containing $cum.start
#' @param dat.type Either "data" or "seg"
#' @importFrom GenomicRanges seqnames
.addCumPos <- function(dat, ref, dat.type){
m.row.idx <- match(as.character(dat$chrom), as.character(seqnames(ref)))
if(dat.type=='data'){
dat$cpos <- ref[m.row.idx,]$cum.start + dat$pos - 1
dat$chr.stat
} else if(dat.type == 'seg'){
dat$cloc.start <- ref[m.row.idx,]$cum.start + dat$loc.start - 1
dat$cloc.end <- ref[m.row.idx,]$cum.start + dat$loc.end - 1
}
dat$chr.stat <- (m.row.idx %% 2) + 1
return(dat)
}
#### Main Functions ####
########################
#' multiDriftPlot OUTDATED
#'
#' @param seg Seg data
#' @param chr.size.gr NULL
#' @param ref.ds NULL
#' @param alt.ds NULL
#' @importFrom methods as
#' @importFrom GenomicRanges seqnames
#' @importFrom GenomicRanges findOverlaps
#' @importFrom GenomicRanges start
#' @importFrom GenomicRanges end
#' @importFrom S4Vectors queryHits
#' @importFrom S4Vectors subjectHits
#' @importFrom graphics plot
#' @importFrom graphics axis
#' @importFrom graphics rect
#'
multiDriftPlot <- function(seg, chr.size.gr=NULL,
ref.ds=NULL, alt.ds=NULL){
if(is.null(ref.ds)) stop("Requires input of ref.ds (GDSC or CCLE)")
if(is.null(alt.ds)) stop("Requires input of alt.ds (GDSC or CCLE)")
if(is.null(chr.size.gr)) chr.size.gr <- .getChrLength()
chrs <- gsub("^chr", "", as.character(seqnames(chr.size.gr)))
# seg <- seg.sig[-null.idx][[2]]
grl <- as(unlist(seg), "GRangesList")
grl.idx <- setNames(c(.grepNA(paste0(ref.ds, "_.*", alt.ds, "_"), names(grl)),
.grepNA(paste0("RNA_.*", alt.ds, "_"), names(grl)),
.grepNA(paste0("RNA_.*", ref.ds, "_"), names(grl))),
c(paste0(ref.ds, "/", alt.ds),
paste0("RNA/", alt.ds),
paste0("RNA/", ref.ds)))
## Initialize the plotting space
plot(0, type='n', xlim=c(1, max(chr.size.gr$cum.end)), ylim=c(0, 3),
xlab="Chr", yaxt='n', ylab='', axes=FALSE)
axis(side=1, at=chr.size.gr$cum.mid[c(TRUE,FALSE)],
labels = chrs[c(TRUE,FALSE)], tick=FALSE, cex.axis=0.8)
axis(side=1, at=chr.size.gr$cum.mid[c(FALSE,TRUE)], line=0.5,
labels = chrs[c(FALSE, TRUE)], tick=FALSE, cex.axis=0.8)
axis(side=2, at=seq(0.5, 2.5, by=1), labels=names(grl.idx), cex.axis=0.8,
las=1, tick=FALSE, line=-1)
lapply(seq_along(grl.idx), function(i){
## Draw empty chromosomes
rect(xleft = chr.size.gr$cum.start, ybottom = rep(i - 0.9, length(chr.size.gr)),
xright = chr.size.gr$cum.end, ytop = rep(i - 0.1, length(chr.size.gr)),
col = "gray95", border="black", lwd=1)
## populate with sig.diff regions
if(!is.na(grl.idx[i])){
gr0 <- grl[[grl.idx[i]]]
ov <- findOverlaps(gr0, chr.size.gr)
gr0$cum.start <- start(gr0[queryHits(ov)]) + chr.size.gr[subjectHits(ov)]$cum.start
gr0$cum.end <- end(gr0[queryHits(ov)]) + chr.size.gr[subjectHits(ov)]$cum.start
with(gr0, rect(xleft = cum.start, ybottom = rep(i - 0.89, length(gr0)),
xright = cum.end, ytop = rep(i - 0.11, length(gr0)),
col = "#fb6a4a", border = NA))
}
})
}
#' plot.CCLid
#' plot() function for CCLid adjusted DNAcopy segment objects
#'
#' @param obj CCLid object returned from bafDrift()
#' @param sample.size Number of points to plot to allow for quicker plotting
#' @param low.sig.alpha low-confident drift regions alpha value for red
#' @param hi.sig.alpha high-confident drift regions alpha value for red
#' @param add.chr.sep Adds lines to delimit chromosomes
#' @param atype default=sd
#' @param add.points Adds points to the segments
#' @param min.z minimum z to report a drift change
#' @importFrom dplyr %>%
#' @import magrittr
#' @importFrom graphics plot
#' @importFrom graphics abline
#' @importFrom graphics axis
#' @importFrom GenomicRanges seqnames
#' @importFrom stats na.omit
#' @export
plot.CCLid <- function (obj, sample.size=600, low.sig.alpha=0.01,
hi.sig.alpha=0.2, add.chr.sep=TRUE,
atype='sd', add.points=TRUE, min.z=3, ...) {
# low.sig.alpha=0.01
# hi.sig.alpha=0.2
# add.chr.sep=TRUE
# sample.size=600
# add.points=TRUE
# atype='sd'
chroms <- paste0("chr", c(1:22, "X", "Y"))
if(any(grepl("(23)|(24)$", obj$data$chrom))){
obj$data$chrom <- gsub("23$", "X", obj$data$chrom) %>%
gsub("24$", "Y", .)
obj$output$chrom <- gsub("23$", "X", obj$output$chrom) %>%
gsub("24$", "Y", .)
}
if(any(!grepl("^chr", obj$data$chrom))){
obj$data$chrom <- paste0("chr", obj$data$chrom)
obj$output$chrom <- paste0("chr", obj$output$chrom)
}
chr.size.dat <- .getChrLength()
chr.data <- .addCumPos(obj$data, chr.size.dat, dat.type='data')
chr.seg <- .addCumPos(obj$output, chr.size.dat, dat.type='seg')
samples <- colnames(obj$data)[-c(1:2)]
chr.cols <- c('black', 'green')
seg.col <- 'orange'
sig.col <- 'red'
ylim <- switch(atype,
"comp"=c(0,0.5),
c(-1,1))
pos.col <- colnames(obj$data)[2]
## Segment yb chromosome
for(s in samples){
size <- min(nrow(chr.data), sample.size)
sample.idx <- sample(1:nrow(chr.data), size=size, replace=FALSE)
sample.dat <- chr.data[sample.idx, , drop=FALSE]
sample.seg <- chr.seg[which(chr.seg$ID %in% s),]
plot(x = sample.dat$cpos, y=sample.dat[,s], col=chr.cols[sample.dat$chr.stat],
pch=16, ylim=ylim, xlim=c(1, max(chr.size.dat$cum.end)), xlab='', ylab=s,
yaxt='n', xaxt='n', axes=FALSE, cex=0.6)
abline(h = 0)
if(add.chr.sep) abline(v = chr.size.dat$cum.end, lwd=0.5, col='grey')
## Adds chr labels
axis(side=1, at = chr.size.dat$cum.mid[c(TRUE, FALSE)],
labels = gsub("^chr", "", as.character(seqnames(chr.size.dat)))[c(TRUE, FALSE)],
las=1, tick=FALSE, lwd = 0, line=-1, cex.axis=0.6)
axis(side=1, at = chr.size.dat$cum.mid[c(FALSE, TRUE)],
labels = gsub("^chr", "", as.character(seqnames(chr.size.dat)))[c(FALSE, TRUE)],
las=1, tick=FALSE, lwd = 0, line=0, cex.axis=0.6)
axis(side = 2, at=seq(-1, 1, by=0.5), labels=seq(-1, 1, by=0.5),
las=1, cex.axis=0.8)
## Adds Signiifcant deviated regions
if(any(na.omit(sample.seg$t) > 0)){
sig.chr.seg <- sample.seg[which(sample.seg$t > 0),]
sig.chr.seg$alpha <- 0
sig.chr.seg$alpha[sig.chr.seg$t < min.z] <- low.sig.alpha
sig.chr.seg$alpha[sig.chr.seg$t >= min.z] <- hi.sig.alpha
with(sig.chr.seg, rect(xleft = cloc.start, ybottom = ylim[1],
xright = cloc.end, ytop = ylim[2],
border=NA, col = scales::alpha(sig.col, alpha)))
}
## Adds Segmean and SD intervals
with(sample.seg, rect(xleft = cloc.start, ybottom = seg.mean - seg.sd,
xright = cloc.end, ytop = seg.mean + seg.sd,
border=NA, col = scales::alpha(seg.col, 0.3)))
with(sample.seg, segments(x0 = cloc.start, y0 = seg.mean,
x1 = cloc.end, y1 = seg.mean,
lwd=3, col=seg.col))
}
}
#' #' plot.SimpleCCLid - OUTDATED
#' #'
#' #' @param obj CCLid object
#' #' @param add.chr.sep Adds chromosome separator
#' #' @param min.z Minimum Z for drift changes
#' plot.SimpleCCLid <- function(obj, add.chr.sep=TRUE, min.z=2){
#' chroms <- paste0("chr", c(1:22, "X", "Y"))
#' if(any(grepl("(23)|(24)$", obj$data$chrom))){
#' obj$data$chrom <- gsub("23$", "X", obj$data$chrom) %>%
#' gsub("24$", "Y", .)
#' obj$output$chrom <- gsub("23$", "X", obj$output$chrom) %>%
#' gsub("24$", "Y", .)
#' }
#' if(any(!grepl("^chr", obj$data$chrom))){
#' obj$data$chrom <- paste0("chr", obj$data$chrom)
#' obj$output$chrom <- paste0("chr", obj$output$chrom)
#' }
#'
#' chr.size.dat <- .getChrLength()
#' chr.seg <- .addCumPos(dat=obj$output, ref=chr.size.dat, dat.type='seg') #CCLid:::
#'
#' seg.col <- c('gray80', 'gray60')
#' sig.col <- 'red'
#'
#' seg.spl <- split(chr.seg, chr.seg$ID)
#'
#' ### Viz begins!
#' ## Create blank canvas to map the brilliance!
#' plot(x = 0, type='n', xlab='', ylab='',
#' pch=16, ylim=c(0, length(seg.spl)), xlim=c(1, max(chr.size.dat$cum.end)),
#' yaxt='n', xaxt='n', axes=FALSE, cex=0.6)
#' abline(h = 0)
#' if(add.chr.sep) abline(v = chr.size.dat$cum.end, lwd=0.5, col='grey')
#'
#' ## Adds chr labels
#' axis(side=1, at = chr.size.dat$cum.mid[c(TRUE, FALSE)],
#' labels = gsub("^chr", "", as.character(seqnames(chr.size.dat)))[c(TRUE, FALSE)],
#' las=1, tick=FALSE, lwd = 0, line=-1, cex.axis=0.6)
#' axis(side=1, at = chr.size.dat$cum.mid[c(FALSE, TRUE)],
#' labels = gsub("^chr", "", as.character(seqnames(chr.size.dat)))[c(FALSE, TRUE)],
#' las=1, tick=FALSE, lwd = 0, line=0, cex.axis=0.6)
#' axis(side = 2, at=(c(1:length(seg.spl)) - 0.5), labels=names(seg.spl),
#' las=1, cex.axis=0.8)
#'
#' for(s in names(seg.spl)){
#' ssid <- seg.spl[[s]]
#' ssid[is.na(ssid)] <- 0
#'
#' ssid$max.t <- ssid$t
#' if(any(ssid$max.t > 5)) ssid[ssid$max.t > 5,]$max.t <- 5
#' y.idx <- match(s, names(seg.spl))
#'
#' ## Add rectangles for the multiplots
#' rect(xleft = ssid$cloc.start, ybottom = (y.idx - 0.8),
#' xright = ssid$cloc.end, ytop = (y.idx - 0.5 + ssid$max.t/10),
#' col=seg.col[ssid$chr.stat], border = NA)
#' if(any(ssid$t >= min.z)){
#' ss.sig <- ssid[which(ssid$t >= min.z),]
#' rect(xleft = ss.sig$cloc.start, ybottom = (y.idx - 0.8),
#' xright = ss.sig$cloc.end, ytop = (y.idx - 0.5 + ss.sig$max.t/10),
#' col=sig.col, border = NA)
#' }
#' }
#'
#' }
#' plot.multiObj OUTDATED
#'
#' @param OBJ A list of CCLid Objects
#' @param add.chr.sep Adds chromosome separator
#' @param min.z Z threshold to report drift change
#'
#' @export
# plot.multiObj <- function(OBJ, add.chr.sep=TRUE, min.z=2){
# if(class(OBJ) == 'list'){
# par(mfrow=c(length(OBJ), 1), mar=c(3, 8, 3, 2))
#
# if(length(min.z) > 1){
# if(length(min.z) != length(OBJ)) stop("min.z length does not equal input length of CCLid objects")
# lapply(seq_along(OBJ), function(idx){
# plot.SimpleCCLid(OBJ[[idx]], min.z[idx], add.chr.sep=add.chr.sep)
# })
# } else {
# lapply(OBJ, plot.SimpleCCLid, min.z=min.z, add.chr.sep=add.chr.sep)
# }
# } else if(class(OBJ)=='CCLid'){
# par(mar=c(3, 8, 3, 2))
# plot.SimpleCCLid(OBJ, min.z=min.z, add.chr.sep=add.chr.sep)
# }
# }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.