Nothing
plotWindows <- function(seqz.window, m.lty = 1, m.lwd = 3,
m.col = "black", q.bg = "lightblue",
log2.plot = FALSE, n.min = 1, xlim, ylim,
add = FALSE, ...) {
if (log2.plot) {
seqz.window[, c(3, 4, 5)] <- log2(seqz.window[, c(3, 4, 5)])
}
if (!add) {
if (missing(xlim))
xlim <- c(seqz.window$start[1], seqz.window$end[nrow(seqz.window)])
if (missing(ylim))
ylim <- c(min(seqz.window$q0, na.rm = TRUE),
max(seqz.window$q1, na.rm = TRUE))
plot(xlim, ylim, type = "n", ...)
}
seqz.window <- seqz.window[seqz.window$N >= n.min, ]
rect(xleft = seqz.window$start, ybottom = seqz.window$q0,
xright = seqz.window$end, ytop = seqz.window$q1,
col = q.bg, border = NA)
segments(y0 = seqz.window$mean, x0 = seqz.window$start,
x1 = seqz.window$end, lty = m.lty, lwd = m.lwd, col = m.col)
}
gc.plot <- function(gc_list, range.gc = NULL, range.depth = NULL, ...) {
n <- gc_list$n
n[n == 0] <- NA
gc <- gc_list$gc
depth <- gc_list$depth
if (length(range.gc) == 2) {
range.gc <- sort(range.gc)
gc.select <- gc <= range.gc[2] & gc >= range.gc[1]
gc <- gc[gc.select]
n <- n[gc.select, ]
}
if (length(range.depth) == 2) {
range.depth <- sort(range.depth)
depth.select <- depth <= range.depth[2] & depth >= range.depth[1]
depth <- depth[depth.select]
n <- n[, depth.select]
}
colorgram(x = gc, y = depth, z = n, ...)
}
gc.summary.plot <- function(gc_list, mean.col = 1, median.col = 2,
scale.subset = 1.5, ...){
mengc <- mean_gc(gc_list)
medgc <- median_gc(gc_list)
max_depth <- round(max(c(mengc, medgc)) * scale.subset, 0)
gc.plot(gc_list, range.depth = c(0, max_depth), ...)
lines(x = gc_list$gc, y = mengc, lwd = 3, col = mean.col)
lines(x = gc_list$gc, y = medgc, lwd = 3, col = median.col)
legend("topright", c("Mean depth", "Median depth"),
col = c(mean.col, median.col), bg = "white", lty = 1, lwd = 3)
}
cp.plot <- function (cp.table, xlab = "Ploidy", ylab = "Cellularity",
zlab = "Scaled rank LPP",
colFn = colorRampPalette(c("white", "lightblue")), ...) {
z <- matrix(rank(cp.table$lpp), nrow = nrow(cp.table$lpp)) /
length(cp.table$lpp)
map <- makecmap(c(0, 1), colFn = colFn, include.lowest = TRUE)
colorgram(x = cp.table$ploidy, y = cp.table$cellularity, z = z,
map = map, las = 1, xlab = xlab, ylab = ylab, zlab = zlab, ...)
}
cp.plot.contours <- function(cp.table, likThresh = c(0.95),
alternative = TRUE, col = palette(), legend.pos = "bottomright",
pch = 18, alt.pch = 3, ...) {
znormsort <- sort(cp.table$lpp, decreasing = TRUE)
znormcumLik <- cumsum(znormsort)
n <- sapply(likThresh, function(x) sum(znormcumLik < x) + 1)
LikThresh <- znormsort[n]
names(LikThresh) <- paste0(likThresh * 100, "%")
contour(x = cp.table$ploidy, y = cp.table$cellularity, z = cp.table$lpp,
levels = znormsort[n], col = col, drawlabels = FALSE,
xlab = "Ploidy", ylab = "Cellularity", ...)
max.xy <- which(cp.table$lpp == max(cp.table$lpp), arr.ind = TRUE)
points(x = cp.table$ploidy[max.xy[, "row"]],
y = cp.table$cellularity[max.xy[, "col"]], pch = pch)
if (alternative == TRUE){
alt.sol <- alternative.cp.solutions(cp.table)
alt.sol <- alt.sol[-1, ]
points(x = alt.sol$ploidy, y = alt.sol$cellularity, pch = alt.pch)
}
if (!is.na(legend.pos)) {
if (alternative == FALSE) {
legend(legend.pos, legend = c(paste("C.R.", names(LikThresh),
sep = " "), "Point estimate"),
col = c(col[1:length(LikThresh)], "black"),
lty = c(rep(1, length(LikThresh)), NA),
pch = c(rep(NA, length(LikThresh)), pch),
border = NA, bty = "n")
} else {
legend(legend.pos, legend = c(paste("C.R.",
names(LikThresh), sep = " "),
"Point estimate", "Alternative solutions"),
col = c(col[1:length(LikThresh)], "black", "black"),
lty = c(rep(1, length(LikThresh)), NA, NA),
pch = c(rep(NA, length(LikThresh)), pch, alt.pch),
border = NA, bty = "n")
}
}
invisible(LikThresh)
}
chromosome.view <- function(baf.windows, ratio.windows, mut.tab = NULL,
segments = NULL, min.N.baf = 1, min.N.ratio = 1e4, main = "",
vlines = FALSE, legend.inset = c(-20 * strwidth("a", units = "figure"), 0),
CNn = 2, cellularity = NULL, ploidy = NULL,
avg.depth.ratio = NULL, model.lwd = 1, model.lty = "24", model.col = 1,
x.chr.space = 10) {
if (is.null(segments)) {
data.model <- NULL
} else {
if ("CNt" %in% colnames(segments)) {
if (length(c(cellularity, ploidy, avg.depth.ratio)) != 3) {
data.model <- NULL
} else {
data.model <- list()
CNt.max <- max(segments$CNt, na.rm = TRUE) + 1
CNt.min <- 0
baf_types <- baf.types.matrix(CNt.min = CNt.min,
CNt.max = CNt.max, CNn = 2)
data.model$baf <- baf.model.points(cellularity = cellularity,
ploidy = ploidy, baf_types = baf_types,
avg.depth.ratio = avg.depth.ratio)
data.model$baf <- data.frame(CNt = baf_types$CNt,
A = baf_types$CNt - baf_types$B, B = baf_types$B,
data.model$baf)
mufreq_types <- mufreq.types.matrix(CNt.min = CNt.min,
CNt.max = CNt.max, CNn = CNn)
data.model$muf <- cbind(mufreq_types,
mufreq.model.points(cellularity = cellularity,
ploidy = ploidy, mufreq_types = mufreq_types,
avg.depth.ratio = avg.depth.ratio))
}
} else {
data.model <- NULL
}
}
if (is.null(mut.tab)) {
par(mar = c(0, 4, 0, 3), oma = c(5, 0, 4, 0),
mfcol = c(2, 1), xaxt = "n")
min.x <- min(c(min(baf.windows$start), min(ratio.windows$start)))
max.x <- max(c(max(baf.windows$end), max(ratio.windows$end)))
xlim <- c(min.x, max.x)
} else {
min.x <- min(c(min(baf.windows$start), min(ratio.windows$start),
min(mut.tab$position)))
max.x <- max(c(max(baf.windows$end), max(ratio.windows$end),
max(mut.tab$position)))
xlim <- c(min.x, max.x)
par(mar = c(0, 4, 0, 10), oma = c(5, 0, 4, 0),
mfcol = c(3, 1), xaxt = "n", xpd = TRUE)
mutation.colors <- c(
"A>C" = rgb(red = 0, green = 178, blue = 238,
alpha = 120, maxColorValue = 255),
"T>G" = rgb(red = 0, green = 178, blue = 238,
alpha = 120, maxColorValue = 255),
"A>G" = rgb(red = 255, green = 64, blue = 64,
alpha = 120, maxColorValue = 255),
"T>C" = rgb(red = 255, green = 64, blue = 64,
alpha = 120, maxColorValue = 255),
"A>T" = rgb(red = 34, green = 139, blue = 34,
alpha = 120, maxColorValue = 255),
"T>A" = rgb(red = 34, green = 139, blue = 34,
alpha = 120, maxColorValue = 255),
"C>A" = rgb(red = 139, green = 90, blue = 0,
alpha = 120, maxColorValue = 255),
"G>T" = rgb(red = 139, green = 90, blue = 0,
alpha = 120, maxColorValue = 255),
"C>G" = rgb(red = 127, green = 0, blue = 255,
alpha = 120, maxColorValue = 255),
"G>C" = rgb(red = 127, green = 0, blue = 255,
alpha = 120, maxColorValue = 255),
"C>T" = rgb(red = 255, green = 215, blue = 0,
alpha = 120, maxColorValue = 255),
"G>A" = rgb(red = 255, green = 215, blue = 0,
alpha = 120, maxColorValue = 255))
plot(x = mut.tab$position, y = mut.tab$F,
ylab = "Mutant allele frequency", las = 1, pch = 19,
col = c(mutation.colors, "NA" = NA)[as.character(mut.tab$mutation)],
ylim = c(min(mut.tab$F, na.rm = TRUE), 1), xlim = xlim)
unique.colors <- unique(mutation.colors)
labels <- sapply(unique.colors, function(a) {
paste(names(mutation.colors)[mutation.colors == a],
collapse = ", ")
})
legend(y = "center", x = "right", legend = labels,
inset = legend.inset, pch = 19, col = unique.colors,
pt.bg = unique.colors, border = NA, bty = "n")
if (!is.null(segments)){
if (vlines) {
abline(v = segments$end.pos, lwd = 1, lty = 2)
}
if (!is.null(data.model)) {
for (i in 1:nrow(segments)) {
segments(x0 = segments$start.pos[i],
x1 = segments$end.pos[i],
y0 = unique(data.model$muf$mufreqs[
data.model$muf$CNt == segments$CNt[i]]),
lwd = model.lwd, lty = model.lty, col = model.col)
}
}
}
}
if (!is.null(segments)){
plot(ylab = "B allele frequency", type = "n",
x = xlim, y = c(0, 0.5), las = 1)
plotWindows(baf.windows, ylab = "B allele frequency",
xlim = xlim, ylim = c(0, 0.5), las = 1,
n.min = min.N.baf, add = TRUE)
if (vlines) {
abline(v = segments$end.pos, lwd = 1, lty = 2)
}
segments(x0 = segments$start.pos, y0 = segments$Bf,
x1 = segments$end.pos, y1 = segments$Bf, col = "red", lwd = 3)
if (!is.null(data.model)) {
for (i in 1:nrow(segments)) {
segments(x0 = segments$start.pos[i], x1 = segments$end.pos[i],
y0 = unique(data.model$baf$BAF[
data.model$baf$CNt == segments$CNt[i]]),
lwd = model.lwd, lty = model.lty, col = model.col)
}
}
} else {
plotWindows(baf.windows, ylab = "B allele frequency",
xlim = xlim, ylim = c(0, 0.5), las = 1, n.min = min.N.baf)
}
plotWindows(ratio.windows, ylab = "Depth ratio",
las = 1, n.min = min.N.ratio, ylim = c(0, 2.5))
if (!is.null(segments)){
if (vlines) {
abline(v = segments$end.pos, lwd = 1, lty = 2)
}
segments(x0 = segments$start.pos, y0 = segments$depth.ratio,
x1 = segments$end.pos, y1 = segments$depth.ratio,
col = "red", lwd = 3)
if (!is.null(data.model)) {
ratios.theoric <- unique(data.model$muf[, c("CNt", "depth.ratio")])
segments(x0 = rep(min(segments$start.pos, na.rm = TRUE),
times = nrow(ratios.theoric)),
x1 = rep(max(segments$end.pos, na.rm = TRUE),
times = nrow(ratios.theoric)),
y0 = ratios.theoric$depth.ratio, lwd = model.lwd,
lty = model.lty, col = model.col)
axis(labels = as.character(ratios.theoric$CNt), side = 4,
line = 0, las = 1, at = ratios.theoric$depth.ratio)
mtext(text = "Copy number", side = 4, line = 2,
cex = par("cex.lab") * par("cex"))
}
}
par(xaxt = "s")
axis(labels = as.character(round(seq(xlim[1] / 1e6, xlim[2] / 1e6,
by = x.chr.space), 0)),
side = 1, line = 0, at = seq(xlim[1], xlim[2], by = 1e6 * x.chr.space),
outer = FALSE, cex = par("cex.axis") * par("cex"))
mtext("Position (Mb)", side = 1, line = 3, outer = FALSE,
cex = par("cex.lab") * par("cex"))
mtext(main, 3, outer = TRUE, cex = par("cex.main") * par("cex"), line = 2)
}
genome.view <- function(seg.cn, info.type = "AB", ...) {
chr.order <- unique(seg.cn$chromosome)
seg.list <- split(x = seg.cn[,
c("chromosome", "start.pos", "end.pos", "A", "B", "CNt")],
f = seg.cn$chromosome)
seg.list <- seg.list[order(order(chr.order))]
seg.max <- lapply(X = seg.list, FUN = function(x) x[nrow(x), "end.pos" ])
seg.pos <- lapply(seg.list, "[", TRUE, c("start.pos", "end.pos"))
seg.max <- cumsum(as.numeric(do.call(rbind, seg.max)))
chr.offset <- 0
for (i in 1:length(seg.pos)){
seg.pos[[i]] <- seg.pos[[i]] + chr.offset
colnames(seg.pos[[i]]) <- c("abs.start", "abs.end")
chr.offset <- seg.max[i]
}
seg.max <- sapply(X = seg.pos, FUN = function(x) x[nrow(x), "abs.end" ])
abs.list <- mapply(cbind, seg.list, seg.pos, SIMPLIFY = FALSE)
abs.segments <- do.call(rbind, abs.list)
if (info.type == "AB") {
na_As <- is.na(abs.segments$A)
max_A <- max(abs.segments$A, na.rm = TRUE)
abs.segments$A[na_As] <- abs.segments$CNt[na_As]
plot(x = c(min(abs.segments$abs.start), max(abs.segments$abs.end)),
y = c(-0.1, (max_A + 0.1)), type = "n",
ylab = "Copy number", xlab = "Position (Mb)",
xaxt = "n", yaxt = "n", xaxs = "i", ...)
axis(labels = 0:max_A, at = 0:max_A, side = 2, line = 0, las = 1)
segments(x0 = abs.segments$abs.start, x1 = abs.segments$abs.end,
y0 = (abs.segments$B - 0.1), y1 = (abs.segments$B - 0.1),
col = "blue", lwd = 5, lend = 1)
segments(x0 = abs.segments$abs.start, x1 = abs.segments$abs.end,
y0 = (abs.segments$A + 0.1), y1 = (abs.segments$A + 0.1),
col = "red", lwd = 5, lend = 1)
} else {
min_CNt <- min(abs.segments$CNt, na.rm = TRUE)
max_CNt <- max(abs.segments$CNt, na.rm = TRUE)
plot(x = c(min(abs.segments$abs.start), max(abs.segments$abs.end)),
y = c(min_CNt, max_CNt), type = "n",
ylab = "Copy number", xlab = "Position (Mb)",
xaxt = "n", yaxt = "n", xaxs = "i", ...)
axis(labels = min_CNt:max_CNt,
at = min_CNt:max_CNt,
side = 2, line = 0, las = 1)
segments(x0 = abs.segments$abs.start, x1 = abs.segments$abs.end,
y0 = abs.segments$CNt, y1 = abs.segments$CNt, col = "red",
lwd = 5, lend = 1)
}
abline(v = c(0, seg.max), lty = 3)
for (i in 1:length(abs.list)){
max.pos <- nrow(abs.list[[i]])
mtext(chr.order[i], side = 3, line = 0,
at = sum(abs.list[[i]]$abs.start[1],
abs.list[[i]]$abs.end[max.pos]) / 2)
}
axis(labels = as.character(round(seq(abs.list[[1]]$start.pos[1] / 1e6,
abs.list[[1]]$end.pos[nrow(abs.list[[1]])] / 1e6, by = 50), 0)),
at = seq(abs.list[[1]]$abs.start[1],
abs.list[[1]]$abs.end[nrow(abs.list[[1]])], by = 5e7),
outer = FALSE, cex = par("cex.axis") * par("cex"), side = 1,
line = 1)
}
plotRawGenome <- function(sequenza.extract, cellularity,
ploidy, CNt.max = 7, main = "", mirror.BAF = TRUE, ...){
max.end <- sapply(sequenza.extract$ratio, FUN = function(x) {
max(x$end, na.rm = TRUE)
})
max.end <- c(0, cumsum(as.numeric(max.end)))
chrs <- names(sequenza.extract$ratio)
coords.names <- (max.end + c(diff(max.end) / 2, 0))[1:length(chrs)]
new.coords <- function(win.list, max.end){
lapply(1:length(win.list), FUN = function(x) {
y <- win.list[[x]]
y$start <- y$start + max.end[x]
y$end <- y$end + max.end[x]
y
}
)}
new.coords.segs <- function(segs, max.end){
lapply(1:length(segs), FUN = function(x) {
y <- segs[[x]]
y$start.pos <- y$start.pos + max.end[x]
y$end.pos <- y$end.pos + max.end[x]
y
}
)}
ratio.new <- new.coords(sequenza.extract$ratio, max.end)
BAF.new <- new.coords(sequenza.extract$BAF, max.end)
segs.new <- do.call(rbind,
new.coords.segs(sequenza.extract$segments, max.end))
avg.depth.ratio <- 1
par(mar = c(1, 4, 0, 3), oma = c(5, 0, 4, 0), mfcol = c(2, 1), ...)
if (mirror.BAF) {
AAF.new <- lapply(BAF.new, function(x) {
x[, 3:5] <- 1 - x[, 3:5]
x
})
plot(x = c(min(max.end), max(max.end)), y = c(0, 1), main = main,
xlab = NA, ylab = "Allele frequency", type = "n", las = 1,
xaxs = "i", yaxs = "i", xaxt = "n" )
plotWindows(seqz.window = do.call(rbind, AAF.new),
q.bg = "lightblue", m.col = "black", add = T)
segments(x0 = segs.new$start.pos, x1 = segs.new$end.pos,
y0 = 1 - (segs.new$Bf), y1 = 1 - (segs.new$Bf),
col = "red", lwd = 2, lend = 1)
} else {
plot(x = c(min(max.end), max(max.end)), y = c(0, 0.5), main = main,
xlab = NA, ylab = "B allele frequency", type = "n", las = 1,
xaxs = "i", yaxs = "i", xaxt = "n" )
}
plotWindows(seqz.window = do.call(rbind, BAF.new), q.bg = "lightblue",
m.col = "black", add = T)
segments(x0 = segs.new$start.pos, x1 = segs.new$end.pos,
y0 = segs.new$Bf, y1 = segs.new$Bf, col = "red", lwd = 2, lend = 1)
abline(v = max.end, lty = 1)
plot(x = c(min(max.end), max(max.end)), y = c(0, 2.5), main = "",
xlab = NA, ylab = "Depth ratio", type = "n", las = 1, xaxs = "i",
yaxs = "i", xaxt = "n")
plotWindows(seqz.window = do.call(rbind, ratio.new), q.bg = "lightblue",
m.col = "black", add = T)
segments(x0 = segs.new$start.pos, x1 = segs.new$end.pos,
y0 = (segs.new$depth.ratio), y1 = (segs.new$depth.ratio),
col = "red", lwd = 2, lend = 1)
if (!missing(ploidy) & !missing(cellularity)){
types <- baf.types.matrix(CNt.min = 0, CNt.max = CNt.max, CNn = 2)
depth.ratios <- baf.model.points(cellularity = cellularity,
ploidy = ploidy, avg.depth.ratio = avg.depth.ratio,
baf_types = types)[, "depth.ratio"]
depth.ratios <- unique(data.frame(CNt = types$CNt,
ratio = depth.ratios))
abline(h = depth.ratios$ratio, lty = 2)
axis(labels = as.character(depth.ratios$CNt), side = 4,
line = 0, las = 1, at = depth.ratios$ratio)
mtext(text = "Copy number", side = 4, line = 2,
cex = par("cex.lab") * par("cex"))
}
abline(v = max.end, lty = 1)
axis(labels = chrs, at = coords.names, side = 1, cex.axis = 1)
}
baf.model.view <- function(cellularity, ploidy, segs,
BAF.space = seq(0.001, 0.5, 0.005), ratio.space = seq(0.01, 2.5, 0.05),
avg.depth.ratio = 1, CNt.max = 7, segment.filter = 3e6, col = "black") {
s.b <- mean(segs$sd.BAF, na.rm = TRUE)
s.r <- mean(segs$sd.ratio, na.rm = TRUE)
l.s <- segs$end.pos - segs$start.pos
s.big <- l.s >= segment.filter
test.values <- expand.grid(Bf = BAF.space, ratio = ratio.space,
KEEP.OUT.ATTRS = FALSE)
both.space <- baf.bayes(Bf = test.values$Bf, CNt.max = CNt.max,
CNt.min = 0, depth.ratio = test.values$ratio,
cellularity = cellularity, ploidy = ploidy,
avg.depth.ratio = avg.depth.ratio, sd.Bf = s.b, weight.Bf = 10,
sd.ratio = s.r, weight.ratio = 10, ratio.priority = F, CNn = 2)
both.space <- as.data.frame(both.space)
z <- tapply(both.space$LPP, list(test.values$Bf, test.values$ratio), mean)
x <- as.numeric(rownames(z))
y <- as.numeric(colnames(z))
t <- baf.types.matrix(CNt.min = 0, CNt.max = CNt.max, CNn = 2)
mpts <- cbind(t, baf.model.points(cellularity = cellularity,
ploidy = ploidy, baf_types = t, avg.depth.ratio = avg.depth.ratio)
)
mpts <- unique(mpts[, c("CNt", "depth.ratio")])
par(mar = c(5.1, 4.1, 4.1, 4.1))
rev.heat <- function(...){rev(heat.colors(...))}
suppressWarnings(colorgram(x, y, z, key = NA, nz = 1000,
xlab = "B allele frequency", ylab = "Depth ratio",
main = paste("cellularity:", cellularity, "ploidy:", ploidy,
"sd.BAF:", round(s.b, 2), sep = " "),
map = makecmap(z, breaks = unique(quantile(z, seq(0.25, 1, 0.0001))),
right = TRUE, n = 1000, colFn = rev.heat), outlier = "white",
las = 1, xlim = c(0, 0.5)))
axis(side = 4, at = mpts$depth.ratio, labels = mpts$CNt, las = 1)
mtext(text = "Copy number", side = 4, line = 2)
segs$col <- col
points(x = segs$Bf[s.big], y = segs$depth.ratio[s.big],
pch = 1, cex = 1, col = segs$col[s.big])
points(x = segs$Bf[!s.big], y = segs$depth.ratio[!s.big], pch = ".",
cex = 1, col = segs$col[!s.big])
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.