addHighlightLegend2 <- function (data, positionsToHighlight, scale)
{
if ((is.null(positionsToHighlight))) {
return()
}
names(positionsToHighlight) = c("chr", "st", "name")
addpts = merge(data, positionsToHighlight, by.x = c("chr",
"st"), by.y = c("chr", "st"))
if (length(addpts[, 1]) == 0) {
return()
}
non.trivial.names <- addpts$name[addpts$name != ""]
if (length(non.trivial.names) == 0) {
return()
}
par(new=TRUE)
#plot.default(x = -10000, y = 1, type = "p", pch = 19, cex = 0.4,
# col = "#00000000", xlim = c(0, 1000), ylim = c(0, 1000),
# axes = FALSE, ann = FALSE)
plot.default(x = -10000, y = 1, type = "p", pch = 19, cex = 0.4,
col = "#00000000", xlim = c(0, 100), ylim = c(0, 1000),
axes = FALSE, ann = FALSE, font=3)
#ypos = rev(seq(0, 900, (900/13)))[1:13]
#step.w = 12
#ypos = rev(seq(0, 900, (800/step.w)))[1:step.w]
step.w = nrow(na.omit(positionsToHighlight))
ypos = rev(seq(20, 900, length.out=step.w))
if (step.w == 1) { ypos = c(100) }
if (step.w == 2) { ypos = rev(seq(100, 200, 100)) }
if (step.w == 3) { ypos = rev(seq(100, 300, 100)) }
if (step.w == 4) { ypos = rev(seq(100, 400, 100)) }
if (step.w == 5) { ypos = rev(seq(100, 500, 100)) }
if (step.w == 6) { ypos = rev(seq(100, 600, 100)) }
if (step.w <= 6) {ypos = ypos - 50}
non.trivial.indices <- (1:length(addpts$name))[addpts$name != ""]
ncol = ceiling(length(non.trivial.names)/step.w)
xpos = 40
offset = 1
nxt <- 1
for (n in 1:ncol) {
names = non.trivial.names[offset:(offset + 12)]
names = as.character(names[!(is.na(names))])
num = length(names)
for (i in 1:num) {
text(xpos, ypos[i], paste(nxt, ". ", names[i], sep = ""),
cex = 0.5, pos = 4, font=3, col='grey40')
nxt <- nxt + 1
}
xpos = xpos + 250
offset = offset + step.w
}
}
plot.SciClone <- function (sco, outputFile, cnToPlot = c(1, 2, 3, 4), showCopyNumberScatterPlots = TRUE,
highlightSexChrs = TRUE, positionsToHighlight = NULL, highlightsHaveNames = FALSE,
overlayClusters = TRUE, overlayIndividualModels = TRUE, showHistogram = FALSE,
showTitle = TRUE, biggerText = FALSE, highlightsOnHistogram = FALSE,
highlightCnPoints = FALSE, smp=NULL)
{
densityData = sco@densities
vafs.merged = sco@vafs.merged
sampleNames = sco@sampleNames
dimensions = sco@dimensions
#nclust = length(unique(sc@clust$cluster.assignments))
nclust = dim(sc@clust$cluster.lower)[2]
if (max(cnToPlot) > 4 | min(cnToPlot) < 1) {
print("sciClone supports plotting of copy numbers between 1 and 4 at this time")
}
minimumLabelledPeakHeight = 999
onlyLabelHighestPeak = TRUE
add.legend <- FALSE
addpts <- NULL
if (!is.null(positionsToHighlight)) {
names(positionsToHighlight) = c("chr", "st", "name")
addpts = merge(vafs.merged, positionsToHighlight, by.x = c("chr",
"st"), by.y = c("chr", "st"))
if ((dim(addpts)[1] > 0) & (any(addpts$name != ""))) {
if (showCopyNumberScatterPlots & (length(cnToPlot) <
2) & highlightsHaveNames) {
add.legend <- TRUE
}
}
}
if (highlightsHaveNames & add.legend) {
if (!(is.null(positionsToHighlight))) {
if (length(positionsToHighlight) < 3) {
print("ERROR: named plot requires names in the third column of the positionsToHighlight data frame")
return(0)
}
cnToPlot = c(2)
}
else {
print("ERROR: highlightsHaveNames requires a 3-column dataframe of positions and names (chr, pos, name)")
return(0)
}
}
clust = NULL
if (overlayClusters) {
if (is.null(sco@clust[1])) {
print("ERROR: can't overlay clusters when clustering was not done on the input data")
return(0)
}
else {
clust = sco@clust
}
}
num.rows <- length(cnToPlot) + 1
if (add.legend) {
# num.rows <- num.rows + 1
}
textScale = 1
axisPosScale = 1
if (biggerText) {
textScale = 1.4
axisPosScale = 0.9
}
height <- 8.5 * (num.rows/5)
width <- 3.7
spacing = 1
scale = 1
if (num.rows == 2) {
spacing = 1.5
scale = 1.5
}
if (num.rows == 3) {
spacing = 1
scale = 1
}
pdf(file = outputFile, width = width, height = height, bg = "white")
numClusters = 0
if (!is.null(clust)) {
numClusters = max(clust$cluster.assignments)
}
for (d in 1:dimensions) {
name = sampleNames[d]
vafs = getOneSampleVafs(vafs.merged, name, numClusters)
par(mfcol = c(num.rows, 1), mar = c(0.5, 3/spacing, 1,
1.5/spacing), oma = c(3/spacing, 0.5, 4/spacing,
0.5), mgp = c(3, 1, 0))
densities = densityData[[d]]$densities
factors = densityData[[d]]$factors
peakPos = densityData[[d]]$peakPos
peakHeights = densityData[[d]]$peakHeights
maxDepth = densityData[[d]]$maxDepth
maxDensity = densityData[[d]]$maxDensity
scalingFactor = 1/maxDensity
plot.default(x = c(1:10), y = c(1:10), ylim = c(0, 1.1),
xlim = c(0, 100), axes = FALSE, ann = FALSE, col = "#00000000",
xaxs = "i", yaxs = "i")
rect(0, 0, 100, 1.1, col = "#00000011", border = NA)
axis(side = 2, at = c(0, 2), labels = c("", ""), las = 1,
cex.axis = 0.6 * textScale, hadj = 0.6 * textScale,
lwd = 0.5, lwd.ticks = 0.5, tck = -0.01)
colors = c("#1C366099", "#67B32E99", "#F4981999", "#E5242099")
density.curve.width <- 4
for (i in cnToPlot) {
if (!(is.null(densities[[i]])) & (!showHistogram |
(i != 2))) {
lines(densities[[i]]$x, scalingFactor * factors[[i]],
col = colors[i], lwd = density.curve.width)
if (length(peakHeights[[i]]) > 0) {
ppos = c()
if (onlyLabelHighestPeak) {
ppos = which(peakHeights[[i]] == max(peakHeights[[i]]))
}
else {
ppos = which((peakHeights[[i]] == max(peakHeights[[i]])) &
(peakHeights[[i]] > minimumLabelledPeakHeight))
}
if (length(ppos) > 0) {
text(x = peakPos[[i]][ppos], y = (scalingFactor *
peakHeights[[i]][ppos]) + 1.7, labels = signif(peakPos[[i]][ppos],
3), cex = 0.7, srt = 0, col = colors[[i]])
}
}
}
else if (showHistogram & (i == 2)) {
v = vafs[which(vafs$cleancn == 2 & vafs$adequateDepth ==
1), ]
frequencies <- data.frame(x = v$vaf, row.names = NULL,
stringsAsFactors = NULL)
bin.width <- 2.5
num.breaks <- ceiling(100/bin.width) + 1
breaks <- unlist(lapply(0:(num.breaks - 1), function(x) 100 *
x/(num.breaks - 1)))
h <- hist(v$vaf, breaks = breaks, plot = FALSE)
h$density <- h$density/max(h$density)
plot(h, add = TRUE, freq = FALSE, col = "white",
border = "black")
}
}
model.style <- 4
model.style <- 1
model.width <- density.curve.width/2
individual.model.style <- 3
individual.model.width <- density.curve.width/2
if (!(is.null(clust))) {
maxFitDensity <- max(clust$fit.y[d, ])
lines(clust$fit.x, clust$fit.y[d, ]/maxFitDensity,
type = "l", col = "grey50", lty = model.style,
lwd = model.width)
if (overlayIndividualModels == TRUE) {
for (i in 1:numClusters) {
lines(clust$fit.x, clust$individual.fits.y[[i]][d,
]/maxFitDensity, type = "l", col = "grey50",
lty = individual.model.style, lwd = individual.model.width)
}
}
if (highlightsOnHistogram) {
if (!is.null(positionsToHighlight)) {
addpts = merge(vafs, positionsToHighlight,
by.x = c("chr", "st"), by.y = c("chr", "st"))
for (i in 1:length(addpts$vaf)) {
if (addpts$name[i] != "") {
vaf <- addpts$vaf[i]
nearest.indx <- which(unlist(lapply(clust$fit.x,
function(x) abs(x - vaf))) == min(abs(clust$fit.x -
vaf)))[1]
vaf.y <- clust$fit.y[d, nearest.indx]/maxFitDensity
label <- as.character(addpts$name[i])
cex <- 1
text(x = vaf, y = vaf.y, label = "*", cex = cex)
text(x = vaf, y = vaf.y + 0.1, label = label,
cex = cex)
}
}
}
}
}
lcol = colors[cnToPlot]
lty = c(1, 1, 1, 1)
lwd = c(2, 2, 2, 2)
pchs = c(NA, NA, NA, NA)
pt.bgs = lcol
leg = c("1 Copy", "2 Copies", "3 Copies", "4 Copies")
leg = leg[cnToPlot]
if ((length(cnToPlot) == 1) & (cnToPlot[1] == 2)) {
if (showHistogram == FALSE) {
lty = c(1)
lwd = c(2)
pchs = c(NA)
pt.bgs = lcol
}
else {
lcol = "black"
lty = c(0)
lwd = c(0)
pchs = c(22)
pt.bgs = "white"
}
}
if (!(is.null(clust))) {
leg = c(leg, "Model Fit")
lcol = c(lcol, "grey50")
lty = c(lty, model.style)
lwd = c(lwd, model.width)
pt.bgs = c(pt.bgs, "grey50")
pchs = c(pchs, NA)
if (overlayIndividualModels == TRUE) {
leg = c(leg, "Component Fits")
lcol = c(lcol, "grey50")
lty = c(lty, individual.model.style)
lwd = c(lwd, 2)
pt.bgs = c(pt.bgs, "grey50")
pchs = c(pchs, NA)
}
}
legend(x = "topright", lwd = lwd, lty = lty, legend = leg,
col = lcol, bty = "n", cex = ((0.6/scale) * textScale),
y.intersp = 1.25, pch = pchs, pt.bg = pt.bgs)
axis(side = 3, at = c(0, 20, 40, 60, 80, 100), labels = c(0,
20, 40, 60, 80, 100), cex.axis = ((0.6/scale) * textScale),
lwd = 0.5, lwd.ticks = 0.5, padj = (((scale * 3.5) -
3.5) + 1.4) * (1/textScale), tck = -0.05)
#mtext("Variant Allele Frequency", adj = 0.5, padj = -3.1 * (1/textScale), cex = 0.6 * (textScale * 0.8), side = 3)
mtext(paste("Variant Allele Frequency (", smp, ", n=", nclust, ")", sep=""), adj = 0.5, padj = -3.1 * (1/textScale), cex = 0.6 * (textScale * 0.8), side = 3)
mtext("Density (a.u.)", side = 2, cex = 0.6 * (textScale *
0.8), padj = -4.2 * (axisPosScale))
if (showTitle) {
title = ""
if (is.null(sampleNames[d])) {
title = "Clonality Plot"
}
else {
title = paste(sampleNames[d], "Clonality Plot", sep = " ")
}
mtext(title, adj = 0.5, padj = -5 * (1/textScale),
cex = 0.65 * textScale, side = 3)
}
if (showCopyNumberScatterPlots) {
for (i in cnToPlot) {
v = vafs[which(vafs$cleancn == i & vafs$adequateDepth ==
1), ]
drawScatterPlot(v, highlightSexChrs, positionsToHighlight,
colors, i, maxDepth, highlightsHaveNames, overlayClusters,
scale, textScale, axisPosScale, labelOnPlot = FALSE,
highlightCnPoints = highlightCnPoints)
axis(side = 1, at = c(0, 20, 40, 60, 80, 100),
labels = c(0, 20, 40, 60, 80, 100), cex.axis = (0.6/scale) *
textScale, lwd = 0.5, lwd.ticks = 0.5, padj = ((-scale *
5) + 5 - 1.4) * (1/textScale), tck = -0.05)
if (length(cnToPlot) < 2 & highlightsHaveNames) {
addHighlightLegend2(v, positionsToHighlight,scale)
}
else {
if (highlightsHaveNames) {
print("WARNING: highlighted point naming is only supported when plotting only CN2 regions (cnToPlot=c(2))")
print("Instead labeling directly on plot")
}
}
}
}
}
devoff <- dev.off()
}
#library(sciClone)
#files = list.files(pattern='.Rdata')
#for (fn in files) {
# load(fn)
# fig=paste(smp,".pdf",sep="")
# plot1d(sc,fig, highlightSexChrs=FALSE, positionsToHighlight=anno.smp, highlightsHaveNames=TRUE, overlayClusters=TRUE, showTitle=FALSE, cnToPlot=c(2), biggerText=FALSE)
# cat(smp, fn, fig, "\n")
#}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.