R/plots.R
In genome/sciclone: Detect subclones from genomic sequence data
Defines functions
getOneSampleVafs
getClusterColors
sc.plot3d
sc.plot2d
compute.binomial.error.bars
sc.plot2dWithMargins
addHighlightLegend
drawScatterPlot
sc.plot1d
Documented in
addHighlightLegend
compute.binomial.error.bars
drawScatterPlot
getClusterColors
getOneSampleVafs
sc.plot1d
sc.plot2d
sc.plot2dWithMargins
sc.plot3d
#---------------------------------------------------------------------------------
## Create the one dimensional plot with kde and scatter
##
sc.plot1d <- function(sco, outputFile=NULL,
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){
densityData = sco@densities
vafs.merged = sco@vafs.merged
sampleNames = sco@sampleNames
dimensions = sco@dimensions
if(max(cnToPlot) > 4 | min(cnToPlot) < 1){
print("sciClone supports plotting of copy numbers between 1 and 4 at this time")
}
#leaving this peak-labeling code in for now, but turn it off
#until it's improved
minimumLabelledPeakHeight=999
onlyLabelHighestPeak=TRUE
# If any of the vafs are named, assume we will be plotting them and
# will need a legend for them.
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(showCopyNumberScatterPlots & ( length(cnToPlot) < 2 ) & highlightsHaveNames) {
add.legend <- TRUE
}
}
#sanity checks
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;
}
## 3.3 x 7.5 is a good dimensionality for 5 rows. Scale accordingly
## if we have fewer rows.
spacing=1.0
scale=1
if(num.rows == 2){
spacing=1.5
scale=1.5
}
if(num.rows == 3){
spacing=1
scale=1
}
height <- 8.5 * (num.rows/5)
width <- 3.7
if(!is.null(outputFile)) {
pdf(file=outputFile, width=width, height=height, bg="white");
}
numClusters = 0
if(!is.null(clust)) {
numClusters = max(clust$cluster.assignments)
}
##one plot for each sample
for(d in 1:dimensions){
name=sampleNames[d]
##grab only the vafs for this sample:
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.0/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
##draw the density plot
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");
##plot bg color
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 for different copy numbers
colors=c("#1C366099","#67B32E99","#F4981999","#E5242099")
density.curve.width <- 4
for(i in cnToPlot){
if(!(is.null(densities[[i]])) & (!showHistogram | (i!= 2))){
##density lines
lines(densities[[i]]$x, scalingFactor*factors[[i]], col=colors[i], lwd=density.curve.width);
##peak labels
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)) {
## Only show histogram for copy number = 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)
## Rescale the intensity so it has max value y = 1.
h$density <- h$density / max(h$density)
plot(h, add=TRUE, freq=FALSE, col="white", border="black")
}
}
##if we have X/Y vals from the clustering alg, add them
## If we overlay the model, show it in a different style
## style 4 = dot dashed
## style 1 = line
model.style <- 4
model.style <- 1
model.width <- density.curve.width / 2
## Plot the individual models with dotted lines (3) or dashed (2)
individual.model.style <- 3
individual.model.width <- density.curve.width / 2
if(!(is.null(clust))){
maxFitDensity <- max(clust$fit.y[d,])
#points(clust$fit.x, clust$fit.y[d,]*25, type="l",col="grey50")
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] != "") {
## We won't evaluate the probability at this VAF--instead
## just look for the closest point at which we evaluated the
## density
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 + .1, label=label, cex=cex)
}
}
}
}
}
##legend
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("Density (a.u.)",side=2,cex=0.6*(textScale*0.8),padj=-4.2*(axisPosScale));
##add a title to the plot
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);
}
##-----------------------------------------------------
##create the scatterplots of vaf vs copy number
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){
addHighlightLegend(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")
}
}
}
}
}
##close the pdf
if(!is.null(outputFile)) {
devoff <- dev.off();
}
}
##--------------------------------------------------------------------
## draw a scatter plot of vaf vs depth
##
drawScatterPlot <- function(data, highlightSexChrs, positionsToHighlight, colors, cn, maxDepth, highlightsHaveNames, overlayClusters, scale=1, textScale=1, axisPosScale=1, labelOnPlot=FALSE, highlightCnPoints=FALSE){
cex.points = 1/scale
## define plot colors
ptcolor = colors[cn]
circlecolor = substr(colors[cn],1,7) #chop off the alpha value
## define the plot space by plotting offscreen points
plot.default( x=-10000, y=1, log="y", type="p", pch=19, cex=0.4,
col="#00000000", xlim=c(0,100), ylim=c(5,maxDepth*3),
axes=FALSE, ann=FALSE, xaxs="i", yaxs="i");
addPoints <- function(data, color, highlightSexChrs, pch=16, cex=0.75, highlightCnPoints){
outlineCol = rgb(0,0,0,0.1);
data.sex=NULL
data.cn=NULL
if(highlightSexChrs){
##plot sex chromsomes with different shape
data.sex = data[(data$chr == "X" | data$chr == "Y"),]
##store autosomes
data = data[!(data$chr == "X" | data$chr == "Y"),]
}
if(highlightCnPoints){
##highlight CN-derived points
data.cn = data[data$chr == "CN",]
##store the rest
data = data[!(data$chr == "CN"),]
}
##plot points normally
points(data$vaf, data$depth, type="p", pch=pch, cex=cex, col=color);
if(!(is.null(data.sex))){
##plot sex chromsomes with different shape
points(data.sex$vaf, data.sex$depth, type="p", pch=pch+1, cex=cex, col=color);
}
if(!(is.null(data.cn))){
##plot cn events with different shape
points(data.cn$vaf, data.cn$depth, type="p", pch=15, cex=cex+0.5, col="black");
points(data.cn$vaf, data.cn$depth, type="p", pch=15, cex=cex, col="yellow");
}
}
#do we have any points to plot?
if(length(data[,1]) > 0){
##if we have cluster assignments in col 8, color them
if(any(grepl(pattern="^cluster$",names(data))) & overlayClusters & cn==2){
numClusters=max(data[,c("cluster")],na.rm=T)
cols = getClusterColors(numClusters)
##plot cluster zero points in grey
p = data[data$cluster == 0,]
if(length(p[,1]) > 0){
addPoints(p, rgb(0,0,0,0.20), highlightSexChrs, cex=(cex.points*0.6), highlightCnPoints=highlightCnPoints)
}
## then the clustered points
for(i in 1:numClusters){
p = data[data$cluster == i,]
addPoints(p, cols[i], highlightSexChrs, cex=cex.points, highlightCnPoints=highlightCnPoints)
}
} else { ##just use the normal color
addPoints(data[data$chr != "CN",], ptcolor, highlightSexChrs, cex=cex.points, highlightCnPoints=highlightCnPoints)
##highlight CN-derived points
addPoints(data[data$chr == "CN",], "black", highlightSexChrs, pch=15, cex=cex.points+0.5, highlightCnPoints=highlightCnPoints)
addPoints(data[data$chr == "CN",], "yellow", highlightSexChrs, pch=15, cex=cex.points, highlightCnPoints=highlightCnPoints)
}
##TODO - add a legend for point types - highlight vs sex chrs vs autosomes
## add highlighted points selected for by user
if(!(is.null(positionsToHighlight))){
addpts = merge(data, positionsToHighlight, by.x=c("chr","st"), by.y = c("chr","st"))
xs <- c()
ys <- c()
labels <- c()
if(length(addpts[,1]) > 0){
if(highlightsHaveNames){
if(labelOnPlot) {
xs <- addpts$vaf[addpts$name != ""]
ys <- addpts$depth[addpts$name != ""]
labels <- addpts$name[addpts$name != ""]
num.labels <- length(labels)
suppressPackageStartupMessages(library(TeachingDemos))
# By using the min argument, ensure that the annotations
# do not overlap the corresponding symbol. NB: we anticipate
# this code will only be active for the interior points
text(x=addpts$vaf,y=addpts$depth,label=rep("*", length(addpts$vaf)))
if(num.labels > 0){
xs <- spread.labs(xs, mindiff=4, min=xs)
ys <- spread.labs(ys+50, mindiff=4, min=ys)
text(x=xs,y=ys,label=labels)
}
} else {
for(i in 1:length(addpts$vaf)){
if(addpts$name[i] != "") {
# Only label the gene with a number if it has a number
text(addpts$vaf[i],addpts$depth[i],labels=i,cex=0.5)
} else {
# Otherwise, just highlight it with a star
text(addpts$vaf[i],addpts$depth[i],labels="*")
}
}
}
} else {
addPoints(addpts, col="#555555FF", highlightSexChrs, cex=cex.points, highlightCnPoints=highlightCnPoints);
}
}
}
}
## define the axis
axis(side=2,las=1,tck=0,lwd=0,cex.axis=1.2/(scale*2)*textScale, hadj=0.5/scale)
for (i in 2:length(axTicks(2)-1)) {
lines(c(-1,101),c(axTicks(2)[i],axTicks(2)[i]),col="#00000022");
}
## plot the background color
rect(-1, 5, 101, axTicks(2)[length(axTicks(2))]*1.05, col = "#00000011",border=NA);
## add cn circle
cnpos = axTicks(2)[length(axTicks(2))-1]
points(x=c(95), y=cnpos, type="p", pch=19, cex=3/scale, col=circlecolor);
text(c(95),y=cnpos, labels=c(cn), cex=1/scale, col="#FFFFFFFF")
## y axis label
mtext("Tumor Coverage", side=2, cex=0.6*(textScale*0.8), padj=-4.2*(axisPosScale));
}
##--------------------------------------------------------------------
## add a legend for highlighted points with names
addHighlightLegend <- 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() }
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);
# mtext("Genes",side=2,cex=0.5,padj=-4.2);
ypos=rev(seq(0,900,(900/13)))[1:13]
# Don't put empty names on legend
non.trivial.indices <- (1:length(addpts$name))[addpts$name != ""]
ncol=ceiling(length(non.trivial.names)/13)
xpos=0;
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(non.trivial.indices[nxt],". ",names[i],sep=""), cex=0.6, pos=4)
text(xpos, ypos[i], paste(nxt,". ",names[i],sep=""), cex=0.6, pos=4)
nxt <- nxt+1
}
xpos=xpos+250;
offset=offset+13;
}
}
##---------------------------------------------------------------------
## create the two dimensional plot with scatter annotated with
## clustering results and 1D plots along margins, this time using
## ggplot2
sc.plot2dWithMargins <- function(sco, outputFile, positionsToHighlight=NULL, highlightsHaveNames=FALSE, overlayErrorBars=FALSE) {
densityData = sco@densities
vafs.merged = sco@vafs.merged
vafs.1d.merged = sco@vafs.1d
sampleNames = sco@sampleNames
dimensions = sco@dimensions
clust = sco@clust
marginalClust = sco@marginalClust
suppressPackageStartupMessages(library(grid))
suppressPackageStartupMessages(library(ggplot2))
plots.1d.list <- list()
res.1d.max.densities <- list()
xmin <- -5
xmax <- 105
# NB: ymin and ymax below on the right-hand side of the = are
# referring to verbatim names in the data argument. But R CMD check
# doesn't understand this, so define these variables. This is only
# a hack to quiet R CMD check and does not affect the ggplot call.
ymin <- NULL
ymax <- NULL
tmp.file <- tempfile("outputFile.tmp")
pdf(file=tmp.file, width=7.2, height=6, bg="white")
# Create (and store) all possible 1D plots
for(d in 1:dimensions){
# Overlay the histogram of the data on the model fit.
#ylab <- "\nDensity\n"
ylab <- "Density (a.u.)"
title <- ""
# Set max posterior density to max of splinefun
limits <- data.frame(x=c(min(marginalClust[[d]]$fit.x), max(marginalClust[[d]]$fit.x)))
f <- splinefun(marginalClust[[d]]$fit.x, marginalClust[[d]]$fit.y[1,,drop=FALSE])
max.posterior.density <- max(unlist(lapply(seq(from=limits$x[1], to=limits$x[2], by=10^-3), f)))
# xlab <- paste("\n", sampleNames[d], "\n", sep="")
xlab <- paste(sampleNames[d], " VAF", sep="")
numClusters = 0
if(!is.null(vafs.1d.merged[[d]]$cluster)) {
numClusters = max(vafs.1d.merged[[d]]$cluster, na.rm=T)
}
vafs = getOneSampleVafs(vafs.1d.merged[[d]], sampleNames[d], numClusters)
# Only show copy number = 2
v = vafs[which(vafs$cleancn==2 & vafs$adequateDepth==1),];
frequencies <- data.frame(x=v$vaf, row.names=NULL, stringsAsFactors=NULL)
g <- ggplot(data = frequencies, aes(x)) + ggtitle(title) + xlab(xlab) + ylab(ylab)
# g <- g + theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank())
g <- g + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank())
g <- g + theme(plot.margin = unit(c(0,0,0,0), "cm"))
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)))
# This is just to appease R CMD check--it does not effect the
# geom_histogram call below.
..ncount..=NULL;
g <- g + geom_histogram(data = frequencies, mapping=aes(x, y=..ncount..*100), fill="white", colour=gray(0.4), breaks=breaks)
# Need to "print" the graph in order to see its maximum y value
# NB: this is redundant at this point, given that I scale the
# histogram to have height 100 above. i.e., max.density will be 100.
tmp <- print(g)
max.density <- max(tmp[["data"]][[1]]$ymax)
res.1d.max.densities[[d]] <- max.density
hline <- data.frame(x = c(0,100), y=c(-5,-5))
g <- g + geom_line(data = hline, aes(x,y))
#vline <- data.frame(y = c(0,max.density * 1.1), x=c(-5,-5))
vline <- data.frame(y = c(0,100), x=c(-5,-5))
g <- g + geom_line(data = vline, aes(x,y))
scale <- max.density / max.posterior.density
# If we actually performed clustering, show it.
if(numClusters > 0) {
f <- splinefun(marginalClust[[d]]$fit.x, scale*marginalClust[[d]]$fit.y[1,,drop=FALSE])
g <- g + stat_function(data = limits, fun=f, mapping=aes(x))
}
g <- g + coord_cartesian(ylim=c(xmin, max.density*1.1), xlim=c(xmin,xmax))
plots.1d.list[[d]] <- g
}
devoff = dev.off()
unlink(tmp.file)
vplayout <- function(x, y) viewport(layout.pos.row = x, layout.pos.col = y, clip="off")
pdf(file=outputFile, width=7.2, height=6, bg="white")
##create a 2d plot for each pairwise combination of samples
for(d1 in 1:(dimensions-1)){
for(d2 in d1:dimensions){
if(d1==d2){
next
}
xlab <- sampleNames[d1]
xlab <- gsub("\\.", " ", xlab)
#xlab <- paste("\n", xlab, "\n", sep="")
xlab <- paste(xlab, " VAF", sep="")
ylab <- sampleNames[d2]
ylab <- gsub("\\.", " ", ylab)
#ylab <- paste("\n", ylab, "\n", sep="")
ylab <- paste(ylab, " VAF", sep="")
numClusters = 0
if(!is.null(vafs.merged$cluster)) {
numClusters = max(vafs.merged$cluster, na.rm=T)
}
vafs1 = getOneSampleVafs(vafs.merged, sampleNames[d1], numClusters);
vafs2 = getOneSampleVafs(vafs.merged, sampleNames[d2], numClusters);
##get only cn2 points
vafs1 = vafs1[which(vafs1$cleancn==2 & vafs1$adequateDepth==1),]
vafs2 = vafs2[which(vafs2$cleancn==2 & vafs2$adequateDepth==1),]
if(numClusters > 0) {
v = merge(vafs1,vafs2,by.x=c("chr","st","cluster"), by.y=c("chr","st","cluster"),suffixes=c(".1",".2"))
# Remove any outliers--these will have cluster assignment 0
v.outlier <- v[v$cluster == 0,]
v <- v[v$cluster != 0,]
} else {
v = merge(vafs1,vafs2,by.x=c("chr","st"), by.y=c("chr","st"),suffixes=c(".1",".2"))
}
v.no.highlight <- v
if(!(is.null(positionsToHighlight))) {
names(positionsToHighlight)=c("chr","st","name");
chr.start.v <- cbind(v[,"chr"], v[,"st"])
chr.start.highlight <- cbind(positionsToHighlight[,1], positionsToHighlight[,2])
v.no.highlight <- v[!(apply(chr.start.v, 1, paste, collapse="$$") %in% apply(chr.start.highlight, 1, paste, collapse="$$")),]
}
title <- ""
# Plot the points that we will not highlight
frequencies.no.highlight <- data.frame(x=v.no.highlight$vaf.1, y=v.no.highlight$vaf.2, row.names=NULL, stringsAsFactors=NULL)
colvec = c()
if(numClusters > 0) {
clusters <- v.no.highlight$cluster
cols=getClusterColors(numClusters)
colvec = cols[clusters]
clusters <- unlist(lapply(clusters, as.character))
frequencies.no.highlight <- data.frame(x=v.no.highlight$vaf.1, y=v.no.highlight$vaf.2, shape=clusters, colour=clusters, row.names=NULL, stringsAsFactors=FALSE)
# NB: shape and colour below on the right-hand side of the = are
# referring to verbatim names in the data argument. But R CMD check
# doesn't understand this, so define these variables. This is only
# a hack to quiet R CMD check and does not affect the ggplot call.
shape <- NULL
colour <- NULL
g <- ggplot(data = frequencies.no.highlight, aes(x=x, y=y)) + ggtitle(title) + xlab(xlab) + ylab(ylab) + geom_point(data = frequencies.no.highlight, aes(x=x, y=y, shape=shape, colour=colour))
# Plot a legend
g <- g + scale_colour_manual(name = "Clusters", labels = 1:numClusters, values = cols[1:numClusters])
g <- g + scale_shape_manual(name = "Clusters", labels = 1:numClusters, values = 1:numClusters)
} else {
g <- ggplot(data = frequencies.no.highlight, aes(x=x, y=y)) + ggtitle(title) + xlab(xlab) + ylab(ylab) + geom_point(data = frequencies.no.highlight, aes(x=x, y=y))
}
if(overlayErrorBars == TRUE) {
err.bars.1 <- compute.binomial.error.bars(v.no.highlight$var.1, v.no.highlight$depth.1) * 100
err.bars.2 <- compute.binomial.error.bars(v.no.highlight$var.2, v.no.highlight$depth.2) * 100
err.df.x <- data.frame(x=v.no.highlight$vaf.1, y=v.no.highlight$vaf.2, xmin=err.bars.1$lb, xmax=err.bars.1$ub)
err.df.y <- data.frame(x=v.no.highlight$vaf.1, y=v.no.highlight$vaf.2, ymin=err.bars.2$lb, ymax=err.bars.2$ub)
if(numClusters > 0) {
g <- g + geom_errorbarh(data = err.df.x, aes(x=x, y=y, xmin=xmin, xmax=xmax), colour=colvec)
g <- g + geom_errorbar(data = err.df.y, aes(x=x, y=y, ymin=ymin, ymax=ymax), colour=colvec)
} else {
g <- g + geom_errorbarh(data = err.df.x, aes(x=x, y=y, xmin=xmin, xmax=xmax))
g <- g + geom_errorbar(data = err.df.y, aes(x=x, y=y, ymin=ymin, ymax=ymax))
}
}
# Now overlay any points that we will highlight
if(!(is.null(positionsToHighlight))) {
# Merge the data and the positions to highlight by chr (col 1)
# and start (col 2)
addpts = merge(v, positionsToHighlight, by.x=c("chr","st"), by.y = c("chr","st"))
if(dim(addpts)[1] > 0) {
frequencies.highlight <- data.frame(x=addpts$vaf.1, y=addpts$vaf.2, row.names=NULL, stringsAsFactors=NULL)
g <- g + geom_point(data = frequencies.highlight, aes(x=x, y=y), shape="*", size=10, colour="black")
if(overlayErrorBars == TRUE) {
err.bars.1 <- compute.binomial.error.bars(addpts$var.1, addpts$depth.1) * 100
err.bars.2 <- compute.binomial.error.bars(addpts$var.2, addpts$depth.2) * 100
err.df.x <- data.frame(x=addpts$vaf.1, y=addpts$vaf.2, xmin=err.bars.1$lb, xmax=err.bars.1$ub)
g <- g + geom_errorbarh(data = err.df.x, aes(x=x, y=y, xmin=xmin, xmax=xmax), colour="black")
err.df.y <- data.frame(x=addpts$vaf.1, y=addpts$vaf.2, ymin=err.bars.2$lb, ymax=err.bars.2$ub)
g <- g + geom_errorbar(data = err.df.y, aes(x=x, y=y, ymin=ymin, ymax=ymax), colour="black")
}
}
}
# g <- g + theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank())
g <- g + theme_bw() + theme(panel.border = element_blank())
g <- g + theme(plot.margin = unit(c(0,0,0,0), "cm"))
# Put the legend in the top right
g <- g + theme(legend.position = c(1,1), legend.justification=c(1,1))
hline <- data.frame(x = c(0,100), y=c(-5,-5))
g <- g + geom_line(data = hline, aes(x,y))
vline <- data.frame(y = c(0,100), x=c(-5,-5))
g <- g + geom_line(data = vline, aes(x,y))
g <- g + coord_cartesian(xlim=c(xmin, xmax), ylim=c(xmin, xmax))
plot.2d <- g
# Code to override clipping
#plot.2d <- ggplot_gtable(ggplot_build(g))
#plot.2d$layout$clip[plot.2d$layout$name=="panel"] <- "off"
#grid.draw(plot.2d)
plot.1d.1 <- plots.1d.list[[d1]]
plot.1d.2 <- plots.1d.list[[d2]]
grid.newpage()
pushViewport(viewport(layout = grid.layout(2, 2), clip="off"))
text.size <- 10
vp <- vplayout(1,1)
vp11 <- vp
plot.2d <- plot.2d + theme(text = element_text(size = text.size))
print(plot.2d, vp=vp)
if(!(is.null(positionsToHighlight))) {
# Merge the data and the positions to highlight by chr (col 1)
# and start (col 2)
addpts = merge(v, positionsToHighlight, by.x=c("chr","st"), by.y = c("chr","st"))
frequencies.highlight <- data.frame(x=addpts$vaf.1, y=addpts$vaf.2, row.names=NULL, stringsAsFactors=NULL)
# Add the labels
if(dim(addpts)[1] > 0){
if(highlightsHaveNames){
# This code adapted from stackoverflow.com/questions/10536396/using-grconvertx-grconverty-in-ggplot2
# Create a new viewport with clipping disabled so we can
# put text outside the plot
depth <- downViewport('panel.3-4-3-4')
pushViewport(dataViewport(xData=c(0,100), yData=c(0,100), clip='off'))
xs <- list()
ys <- list()
labels <- list()
nxt <- 1
for(i in 1:dim(addpts)[1]) {
if(addpts$vaf.1[i] < 1) {
x <- addpts$vaf.1[i] - 12
y <- addpts$vaf.2[i]
} else if(addpts$vaf.2[i] < 1) {
x <- addpts$vaf.1[i]
y <- addpts$vaf.2[i] - 10
} else {
x <- addpts$vaf.1[i] + 5
y <- addpts$vaf.2[i] + 5
}
label <- as.character(addpts$name[i])
#df <- data.frame(x=x, y=y)
if(label == "") {
# No label to place
} else {
# grid.text(x=x,y=y,label=label,default.units="native", gp=gpar(fontsize=8))
xs[[nxt]] <- x
ys[[nxt]] <- y
labels[[nxt]] <- label
nxt <- nxt + 1
}
}
xs <- unlist(xs)
ys <- unlist(ys)
labels <- unlist(labels)
num.labels <- length(labels)
suppressPackageStartupMessages(library(TeachingDemos))
# By using the min argument, ensure that the annotations
# do not overlap the corresponding symbol. NB: we anticipate
# this code will only be active for the interior points
if(num.labels > 0){
xs <- spread.labs(xs, mindiff=4, min=xs)
ys <- spread.labs(ys, mindiff=4, min=ys)
for(i in 1:num.labels) {
grid.text(x=xs[i],y=ys[i],label=labels[i],default.units="native", gp=gpar(fontsize=8))
}
}
# Move up depth+1 in the tree. NB: +1 because we pushed a
# viewport on to the tree.
upViewport(depth+1)
}
}
}
vp <- vplayout(1,2)
plot.1d.2 <- plot.1d.2 + theme(text = element_text(size = text.size))
plot.1d.2 <- plot.1d.2 + coord_flip(ylim = c(xmin, res.1d.max.densities[[d2]]*1.1), xlim = c(xmin, xmax))
print(plot.1d.2, vp=vp)
vp <- vplayout(2,1)
plot.1d.1 <- plot.1d.1 + theme(text = element_text(size = text.size))
print(plot.1d.1, vp=vp)
}
}
devoff = dev.off()
}
## Compute the lower and upper bound for a "1-std dev" binomial confidence
## interval for a given number of successes and total number trials
compute.binomial.error.bars <- function(successes, total.trials){
suppressPackageStartupMessages(library(MKmisc))
suppressPackageStartupMessages(library(NORMT3))
# Return a "1 std dev" confidence interval
width <- as.double(erf(1/sqrt(2)))
lb <- mapply(function(a,b) binomCI(a, b, conf.level=width, method="jeffreys")$CI[1], successes, total.trials)
ub <- mapply(function(a,b) binomCI(a, b, conf.level=width, method="jeffreys")$CI[2], successes, total.trials)
res <- data.frame(lb=lb, ub=ub)
return(res)
}
##---------------------------------------------------------------------
## create the two dimensional plot with scatter annotated with
## clustering results and 1D plots along margins
##---------------------------------------------------------------------------------
## Create two dimensional plot with scatter annotated with clustering result
##
sc.plot2d <- function(sco, outputFile=NULL, positionsToHighlight=NULL, highlightsHaveNames=FALSE, overlayClusters=TRUE, overlayErrorBars=FALSE, ellipse.metadata = list(), singlePage=FALSE, scale=1, xlim=100, ylim=100, plot.title=NULL, samplesToPlot=NULL, clusterLegend=TRUE, flipSamples=FALSE, xlab=NULL, ylab=NULL, colors=NULL, plotWidth=7.2, plotHeight=6){
vafs.merged = sco@vafs.merged
sampleNames = sco@sampleNames
dimensions = sco@dimensions
if(singlePage){
nplots = ncol(combn(c(1:dimensions),2))
nrow=round(sqrt(nplots))
ncol=ceiling(sqrt(nplots))
if(!is.null(outputFile)){
pdf(outputFile, width=plotWidth*ncol, height=plotHeight*nrow, bg="white")
}
par(mfrow=c(nrow,ncol), mar=c(5.1, 5.1, 4.1, 2.1))
} else {
if(!is.null(outputFile)){
pdf(outputFile, width=plotWidth, height=plotHeight, bg="white")
}
}
numClusters = 0
maxCluster = 0
non.empty.clusters <- c()
# NB: clusters no longer need be numbered contiguously
if(!is.null(vafs.merged$cluster)) {
# Careful! This doesn't work--0 is an outlier cluster.
#numClusters = length(unique(vafs.merged$cluster[!is.na(vafs.merged$cluster)]))
maxCluster = max(vafs.merged$cluster, na.rm=T)
# Note that by starting at 1, we exclude the outlier "cluster" 0.
for(i in 1:maxCluster){
if(!any(vafs.merged$cluster == i)) { next }
non.empty.clusters <- c(non.empty.clusters, i)
}
}
numClusters = length(non.empty.clusters)
suppressPackageStartupMessages(library(plotrix)) # For plotCI among others.
createPlot <- function(d1, d2){
vafs1 = getOneSampleVafs(vafs.merged, sampleNames[d1], numClusters);
vafs2 = getOneSampleVafs(vafs.merged, sampleNames[d2], numClusters);
##get only cn2 points with adequate coverage
vafs1 = vafs1[which(vafs1$cleancn==2 & vafs1$adequateDepth==1),]
vafs2 = vafs2[which(vafs2$cleancn==2 & vafs2$adequateDepth==1),]
if(!is.null(vafs.merged$cluster)) {
v = merge(vafs1,vafs2,by.x=c("chr","st","cluster"), by.y=c("chr","st","cluster"),suffixes=c(".1",".2"))
## Remove any outliers--these will have cluster assignment 0
v.outlier <- v[v$cluster == 0,]
v <- v[v$cluster != 0,]
} else {
v = merge(vafs1,vafs2,by.x=c("chr","st"), by.y=c("chr","st"),suffixes=c(".1",".2"))
}
cols=colors
if(is.null(colors)){
cols = getClusterColors(maxCluster)
}
##sample name
title=paste(sampleNames[d1],"vs",sampleNames[d2])
if(!is.null(plot.title)){
title=plot.title
}
##create the plot
xlabel = paste(sampleNames[d1],"VAF ")
ylabel = paste(sampleNames[d2],"VAF")
if(!(is.null(xlab))){
xlabel = xlab
}
if(!(is.null(ylab))){
ylabel = ylab
}
plot(-100, -100, xlim=c(0,xlim*1.2), ylim=c(0,ylim), main=title,
xlab=xlabel, ylab=ylabel,
bty="n", xaxt="n", yaxt="n", cex.lab=scale, cex.main=scale, cex.axis=scale)
xGridIncrement=20
yGridIncrement=20
if(xlim < 80){
xGridIncrement= 10
}
if(xlim < 50){
xGridIncrement=5
}
if(ylim < 80){
yGridIncrement=10
}
if(ylim < 50){
yGridIncrement=5
}
## vertical grid
abline(v=seq(0,xlim,xGridIncrement),col="grey50", lty=3, lwd=scale)
axis(side=1,at=seq(0,xlim,xGridIncrement),labels=seq(0,xlim,xGridIncrement), cex.axis=scale)
## horizontal grid
numYGridLines=floor(ylim/yGridIncrement)
segments(rep(-10,numYGridLines),seq(0,ylim,yGridIncrement),rep(xlim*1.05,numYGridLines),seq(0,ylim,yGridIncrement), lty=3, lwd=scale, col="grey50")
axis(side=2,at=seq(0,ylim,yGridIncrement),labels=seq(0,ylim,yGridIncrement), cex.axis=scale)
## If we will be highlighting some points, exclude them from
## the general list of points to plot and plot them instead with
## a different symbol/color (a black *)
v.no.highlight <- v
if(!(is.null(positionsToHighlight))) {
names(positionsToHighlight)=c("chr","st","name");
chr.start.v <- cbind(v[,"chr"], v[,"st"])
chr.start.highlight <- cbind(positionsToHighlight[,1], positionsToHighlight[,2])
v.no.highlight <- v[!(apply(chr.start.v, 1, paste, collapse="$$") %in% apply(chr.start.highlight, 1, paste, collapse="$$")),]
}
if(overlayErrorBars == TRUE) {
err.bars.1 <- compute.binomial.error.bars(v.no.highlight$var.1, v.no.highlight$depth.1) * 100
err.bars.2 <- compute.binomial.error.bars(v.no.highlight$var.2, v.no.highlight$depth.2) * 100
}
if(!is.null(vafs.merged$cluster)) {
## handle cluster 0 (outliers)
if(length(v.outlier[,1]) > 0){
points(v.outlier$vaf.1, v.outlier$vaf.2, col=rgb(0,0,0,0.5), pch=".", cex=3*scale)
}
for(i in 1:numClusters){
indices <- v.no.highlight$cluster==non.empty.clusters[i]
if(overlayClusters){
if(dim(v.no.highlight[indices,])[1] > 0) {
if(overlayErrorBars == TRUE) {
plotCI(v.no.highlight[indices,]$vaf.1, v.no.highlight[indices,]$vaf.2, col=cols[non.empty.clusters[i]], pch=i, li=err.bars.1$lb[indices], ui=err.bars.1$ub[indices], add=TRUE, err="x")
plotCI(v.no.highlight[indices,]$vaf.1, v.no.highlight[indices,]$vaf.2, col=cols[non.empty.clusters[i]], pch=i, li=err.bars.2$lb[indices], ui=err.bars.2$ub[indices], add=TRUE, err="y")
} else {
points(v.no.highlight[indices,]$vaf.1, v.no.highlight[indices,]$vaf.2, col=cols[non.empty.clusters[i]], pch=i, cex=scale)
}
}
} else { ##no overlay of clusters
if(dim(v.no.highlight[indices,])[1] > 0) {
if(overlayErrorBars == TRUE) {
plotCI(v.no.highlight[indices,]$vaf.1, v.no.highlight[indices,]$vaf.2, pch=18, li=err.bars.1$lb[indices], ui=err.bars.1$ub[indices], add=TRUE, err="x")
plotCI(v.no.highlight[indices,]$vaf.1, v.no.highlight[indices,]$vaf.2, pch=18, li=err.bars.2$lb[indices], ui=err.bars.2$ub[indices], add=TRUE, err="y")
} else {
points(v.no.highlight[indices,]$vaf.1, v.no.highlight[indices,]$vaf.2, pch=18, cex=scale,col=getClusterColors(1)[1])
}
}
}
}
} else { #no clusters plotted
if(dim(v.no.highlight)[1] > 0) {
if(overlayErrorBars == TRUE) {
plotCI(v.no.highlight$vaf.1, v.no.highlight$vaf.2, pch=18, li=err.bars.1$lb, ui=err.bars.1$ub, add=TRUE, err="x", cex=scale, col=getClusterColors(1)[1])
plotCI(v.no.highlight$vaf.1, v.no.highlight$vaf.2, pch=18, li=err.bars.2$lb, ui=err.bars.2$ub, add=TRUE, err="y", cex=scale, col=getClusterColors(1)[1])
} else {
points(v.no.highlight$vaf.1, v.no.highlight$vaf.2, pch=18, cex=scale, col=getClusterColors(1)[1])
}
}
}
## Now plot the highlighted points so they are overlaid
if(!(is.null(positionsToHighlight))) {
## Merge the data and the positions to highlight by chr (col 1)
## and start (col 2)
addpts = merge(v, positionsToHighlight, by.x=c("chr","st"), by.y = c("chr","st"))
## Plot the highlighted items. NB: we never overlay the
## cluster on them, but expect this will be obvious from context
if(dim(addpts)[1] > 0) {
if(overlayErrorBars == TRUE) {
err.bars.1 <- compute.binomial.error.bars(addpts$var.1, addpts$depth.1) * 100
err.bars.2 <- compute.binomial.error.bars(addpts$var.2, addpts$depth.2) * 100
plotCI(addpts$vaf.1, addpts$vaf.2, pch="*", col="black", cex=2, li=err.bars.1$lb, ui=err.bars.1$ub, add=TRUE, err="x")
plotCI(addpts$vaf.1, addpts$vaf.2, pch="*", col="black", cex=2, li=err.bars.2$lb, ui=err.bars.2$ub, add=TRUE, err="y")
} else {
points(addpts$vaf.1, addpts$vaf.2, pch="*", col="black", cex=2*scale)
}
}
}
if(!is.null(vafs.merged$cluster)) {
for(i in 1:numClusters){
if((!is.null(ellipse.metadata$SEMs.lb)) & (!is.null(ellipse.metadata$SEMs.ub))) {
xc <- ellipse.metadata$SEMs.lb[i,d1] + ((ellipse.metadata$SEMs.ub[i,d1] - ellipse.metadata$SEMs.lb[i,d1])/2)
yc <- ellipse.metadata$SEMs.lb[i,d2] + ((ellipse.metadata$SEMs.ub[i,d2] - ellipse.metadata$SEMs.lb[i,d2])/2)
## ell <- my.ellipse(hlaxa = ((ellipse.metadata$SEMs.ub[i,d1] - ellipse.metadata$SEMs.lb[i,d1])/2), hlaxb = ((ellipse.metadata$SEMs.ub[i,d2] - ellipse.metadata$SEMs.lb[i,d2])/2), xc = xc, yc = yc)
draw.ellipse(xc, yc, a = ((ellipse.metadata$SEMs.ub[i,d1] - ellipse.metadata$SEMs.lb[i,d1])/2), b = ((ellipse.metadata$SEMs.ub[i,d2] - ellipse.metadata$SEMs.lb[i,d2])/2))
}
if((!is.null(ellipse.metadata$std.dev.lb)) & (!is.null(ellipse.metadata$std.dev.ub))) {
xc <- ellipse.metadata$std.dev.lb[i,d1] + ((ellipse.metadata$std.dev.ub[i,d1] - ellipse.metadata$std.dev.lb[i,d1])/2)
yc <- ellipse.metadata$std.dev.lb[i,d2] + ((ellipse.metadata$std.dev.ub[i,d2] - ellipse.metadata$std.dev.lb[i,d2])/2)
## Plot std dev as dashed line.
draw.ellipse(xc, yc, a = ((ellipse.metadata$std.dev.ub[i,d1] - ellipse.metadata$std.dev.lb[i,d1])/2), b = ((ellipse.metadata$std.dev.ub[i,d2] - ellipse.metadata$std.dev.lb[i,d2])/2), lty=2)
}
}
}
##plot(-100, -100, xlim=c(0,100), ylim=c(0,100),main="Clusters")
if(!is.null(vafs.merged$cluster)) {
## Only include in the legend any clusters that are not empty.
if(clusterLegend==TRUE){
legend("topright", legend=non.empty.clusters, col=cols[non.empty.clusters], title="Clusters", pch=non.empty.clusters, cex=scale, bg="white")
}
##legend("topright", legend=1:numClusters, col=cols[1:numClusters], title="Clusters", pch=1:numClusters)
}
## Add annotation for gene names, if requested
if(!(is.null(positionsToHighlight))){
addpts = merge(v, positionsToHighlight, by.x=c("chr","st"), by.y = c("chr","st"))
## write.table(file="genes.txt", unique(addpts$gene_name), row.names=FALSE, col.names=FALSE, quote=FALSE)
if(dim(addpts)[1] > 0){
if(highlightsHaveNames){
xs <- list()
ys <- list()
labels <- list()
nxt <- 1
for(i in 1:dim(addpts)[1]) {
par(xpd=NA)
cex <- 1
if(addpts$vaf.1[i] < 1) {
##text(addpts$vaf.1[i] - 8,addpts$vaf.2[i],labels=addpts$gene_name[i],cex=cex)
x <- addpts$vaf.1[i] - 8
y <- addpts$vaf.2[i]
} else if(addpts$vaf.2[i] < 1) {
##text(addpts$vaf.1[i],addpts$vaf.2[i] - 5,labels=addpts$gene_name[i],cex=cex)
x <- addpts$vaf.1[i]
y <- addpts$vaf.2[i] - 5
} else {
##text(addpts$vaf.1[i] + 4,addpts$vaf.2[i] + 4,labels=addpts$gene_name[i],cex=cex)
x <- addpts$vaf.1[i] + 4
y <- addpts$vaf.2[i] + 4
}
label <- as.character(addpts$name[i])
##df <- data.frame(x=x, y=y)
##grid.text(x=x,y=y,label=label,default.units="native", gp=gpar(fontsize=8))
if(label != "") {
xs[[nxt]] <- x
ys[[nxt]] <- y
labels[[nxt]] <- label
nxt <- nxt + 1
}
}
xs <- unlist(xs)
ys <- unlist(ys)
labels <- unlist(labels)
num.labels <- length(labels)
suppressPackageStartupMessages(library(TeachingDemos))
## By using the min argument, ensure that the annotations
## do not overlap the corresponding symbol. NB: we anticipate
## this code will only be active for the interior points
if(num.labels > 0){
xs <- spread.labs(xs, mindiff=4, min=xs)
ys <- spread.labs(ys, mindiff=4, min=ys)
for(i in 1:num.labels) {
##grid.text(x=xs[i],y=ys[i],label=labels[i],default.units="native", gp=gpar(fontsize=8))
text(x=xs[i], y=ys[i], label=labels[i], cex=cex)
}
}
}
}
} # End add gene annotations
}
##create 2d plots
count=0;
for(d1 in 1:(dimensions-1)){
for(d2 in d1:dimensions){
if(d1==d2){
next
}
if(is.null(samplesToPlot)){ ##use each pairwise combination of samples
if(flipSamples == TRUE){
createPlot(d2,d1)
} else {
createPlot(d1,d2)
}
count = count + 1;
} else { ##just plot the specified samples
if( (sampleNames[d1] %in% samplesToPlot) & (sampleNames[d2] %in% samplesToPlot)){
if(flipSamples == TRUE){
createPlot(d2,d1)
} else {
createPlot(d1,d2)
}
count = count + 1;
}
}
}
}
if(count < 1){
print ("WARNING: no pair of samples matched the samplesToPlot argument. Nothing plotted");
}
if(!is.null(outputFile)){
devoff = dev.off()
}
}
##-------------------------------------------------------------------------------------
## plot three samples in 3d, optionally create a GIF
##
sc.plot3d <- function(sco, samplesToPlot, size=700, outputFile=NULL){
suppressPackageStartupMessages(library(rgl))
##set the size of the window
open3d(windowRect = c(0,0, size, size) )
if(length(samplesToPlot) != 3){
print("ERROR: must provide exactly 3 sample names for 3d plotting")
return(0)
}
a = sco@vafs.merged[,c(paste(samplesToPlot,".vaf",sep=""),"cluster")]
a = a[!is.na(a$cluster),]
a = a[!(a$cluster==0),]
numClusters=max(a$cluster)
cols=getClusterColors(numClusters)
colvec = cols[a$cluster]
plot3d(a[,1], a[,2], a[,3], xlim=c(0,100), ylim=c(0,100),zlim=c(0,100), axes=FALSE,
xlab=samplesToPlot[1], ylab=samplesToPlot[2], zlab=samplesToPlot[3],
type="s", col=colvec)
##add a box
axes3d( edges=c("x--", "y--", "z"),labels=FALSE)
for(i in c("+","-")){
for(j in c("+","-")){
axes3d( edges=paste("x",i,j,sep=""), tick=FALSE, labels=FALSE)
axes3d( edges=paste("y",i,j,sep=""), tick=FALSE, labels=FALSE)
axes3d( edges=paste("z",i,j,sep=""), tick=FALSE, labels=FALSE)
}
}
if(is.null(outputFile)){
play3d(spin3d(axis=c(0,0,1), rpm=10), duration=6)
} else {
##remove trailing .gif, since movie3d adds it
outputFile = sub(".gif$","",outputFile)
movie3d(spin3d(axis=c(0,0,1), rpm=5), duration=12, dir=getwd(), movie=outputFile)
}
}
##-------------------------------------------------------------------------------------
## Get a list of colors to use for the clusters (after 20 they start repeating)
##
getClusterColors <- function(numClusters){
suppressPackageStartupMessages(library(RColorBrewer))
cols=brewer.pal(8,"Dark2")
if(numClusters>8){
cols = c(cols,brewer.pal(9,"Set1"))
}
if(numClusters > 20){ #if this happens, something is probably weird with your data, but we'll support it anyway
cols = rep(cols,ceiling(length(cols)/20))
}
cols = cols[1:numClusters]
##add transparency to colors
for(i in 1:length(cols)){
z = col2rgb(cols[i])
cols[i] = rgb(z[1], z[2], z[3], 150, maxColorValue=255)
}
return(cols)
}
##-------------------------------------------------------------------------------------
## Extract a single sample's data from the merged data frame
##
getOneSampleVafs <- function(vafs.merged, name, numClusters){
common = c("chr","st","adequateDepth")
a = vafs.merged[,common]
prefix = paste("^",name,".",sep="")
cols=grepl(prefix, names(vafs.merged))
header=sub(prefix, "", names(vafs.merged)[cols])
vafs = cbind(a,vafs.merged[,cols])
names(vafs) = c(common,header)
if(numClusters > 0){
vafs = cbind(vafs, vafs.merged$cluster)
names(vafs)[length(vafs)] = "cluster"
}
return(vafs)
}
genome/sciclone documentation built on Oct. 15, 2023, 5:09 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions getOneSampleVafs getClusterColors sc.plot3d sc.plot2d compute.binomial.error.bars sc.plot2dWithMargins addHighlightLegend drawScatterPlot sc.plot1d
Documented in addHighlightLegend compute.binomial.error.bars drawScatterPlot getClusterColors getOneSampleVafs sc.plot1d sc.plot2d sc.plot2dWithMargins sc.plot3d
#---------------------------------------------------------------------------------
## Create the one dimensional plot with kde and scatter
##
sc.plot1d <- function(sco, outputFile=NULL,
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){
densityData = sco@densities
vafs.merged = sco@vafs.merged
sampleNames = sco@sampleNames
dimensions = sco@dimensions
if(max(cnToPlot) > 4 | min(cnToPlot) < 1){
print("sciClone supports plotting of copy numbers between 1 and 4 at this time")
}
#leaving this peak-labeling code in for now, but turn it off
#until it's improved
minimumLabelledPeakHeight=999
onlyLabelHighestPeak=TRUE
# If any of the vafs are named, assume we will be plotting them and
# will need a legend for them.
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(showCopyNumberScatterPlots & ( length(cnToPlot) < 2 ) & highlightsHaveNames) {
add.legend <- TRUE
}
}
#sanity checks
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;
}
## 3.3 x 7.5 is a good dimensionality for 5 rows. Scale accordingly
## if we have fewer rows.
spacing=1.0
scale=1
if(num.rows == 2){
spacing=1.5
scale=1.5
}
if(num.rows == 3){
spacing=1
scale=1
}
height <- 8.5 * (num.rows/5)
width <- 3.7
if(!is.null(outputFile)) {
pdf(file=outputFile, width=width, height=height, bg="white");
}
numClusters = 0
if(!is.null(clust)) {
numClusters = max(clust$cluster.assignments)
}
##one plot for each sample
for(d in 1:dimensions){
name=sampleNames[d]
##grab only the vafs for this sample:
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.0/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
##draw the density plot
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");
##plot bg color
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 for different copy numbers
colors=c("#1C366099","#67B32E99","#F4981999","#E5242099")
density.curve.width <- 4
for(i in cnToPlot){
if(!(is.null(densities[[i]])) & (!showHistogram | (i!= 2))){
##density lines
lines(densities[[i]]$x, scalingFactor*factors[[i]], col=colors[i], lwd=density.curve.width);
##peak labels
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)) {
## Only show histogram for copy number = 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)
## Rescale the intensity so it has max value y = 1.
h$density <- h$density / max(h$density)
plot(h, add=TRUE, freq=FALSE, col="white", border="black")
}
}
##if we have X/Y vals from the clustering alg, add them
## If we overlay the model, show it in a different style
## style 4 = dot dashed
## style 1 = line
model.style <- 4
model.style <- 1
model.width <- density.curve.width / 2
## Plot the individual models with dotted lines (3) or dashed (2)
individual.model.style <- 3
individual.model.width <- density.curve.width / 2
if(!(is.null(clust))){
maxFitDensity <- max(clust$fit.y[d,])
#points(clust$fit.x, clust$fit.y[d,]*25, type="l",col="grey50")
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] != "") {
## We won't evaluate the probability at this VAF--instead
## just look for the closest point at which we evaluated the
## density
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 + .1, label=label, cex=cex)
}
}
}
}
}
##legend
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("Density (a.u.)",side=2,cex=0.6*(textScale*0.8),padj=-4.2*(axisPosScale));
##add a title to the plot
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);
}
##-----------------------------------------------------
##create the scatterplots of vaf vs copy number
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){
addHighlightLegend(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")
}
}
}
}
}
##close the pdf
if(!is.null(outputFile)) {
devoff <- dev.off();
}
}
##--------------------------------------------------------------------
## draw a scatter plot of vaf vs depth
##
drawScatterPlot <- function(data, highlightSexChrs, positionsToHighlight, colors, cn, maxDepth, highlightsHaveNames, overlayClusters, scale=1, textScale=1, axisPosScale=1, labelOnPlot=FALSE, highlightCnPoints=FALSE){
cex.points = 1/scale
## define plot colors
ptcolor = colors[cn]
circlecolor = substr(colors[cn],1,7) #chop off the alpha value
## define the plot space by plotting offscreen points
plot.default( x=-10000, y=1, log="y", type="p", pch=19, cex=0.4,
col="#00000000", xlim=c(0,100), ylim=c(5,maxDepth*3),
axes=FALSE, ann=FALSE, xaxs="i", yaxs="i");
addPoints <- function(data, color, highlightSexChrs, pch=16, cex=0.75, highlightCnPoints){
outlineCol = rgb(0,0,0,0.1);
data.sex=NULL
data.cn=NULL
if(highlightSexChrs){
##plot sex chromsomes with different shape
data.sex = data[(data$chr == "X" | data$chr == "Y"),]
##store autosomes
data = data[!(data$chr == "X" | data$chr == "Y"),]
}
if(highlightCnPoints){
##highlight CN-derived points
data.cn = data[data$chr == "CN",]
##store the rest
data = data[!(data$chr == "CN"),]
}
##plot points normally
points(data$vaf, data$depth, type="p", pch=pch, cex=cex, col=color);
if(!(is.null(data.sex))){
##plot sex chromsomes with different shape
points(data.sex$vaf, data.sex$depth, type="p", pch=pch+1, cex=cex, col=color);
}
if(!(is.null(data.cn))){
##plot cn events with different shape
points(data.cn$vaf, data.cn$depth, type="p", pch=15, cex=cex+0.5, col="black");
points(data.cn$vaf, data.cn$depth, type="p", pch=15, cex=cex, col="yellow");
}
}
#do we have any points to plot?
if(length(data[,1]) > 0){
##if we have cluster assignments in col 8, color them
if(any(grepl(pattern="^cluster$",names(data))) & overlayClusters & cn==2){
numClusters=max(data[,c("cluster")],na.rm=T)
cols = getClusterColors(numClusters)
##plot cluster zero points in grey
p = data[data$cluster == 0,]
if(length(p[,1]) > 0){
addPoints(p, rgb(0,0,0,0.20), highlightSexChrs, cex=(cex.points*0.6), highlightCnPoints=highlightCnPoints)
}
## then the clustered points
for(i in 1:numClusters){
p = data[data$cluster == i,]
addPoints(p, cols[i], highlightSexChrs, cex=cex.points, highlightCnPoints=highlightCnPoints)
}
} else { ##just use the normal color
addPoints(data[data$chr != "CN",], ptcolor, highlightSexChrs, cex=cex.points, highlightCnPoints=highlightCnPoints)
##highlight CN-derived points
addPoints(data[data$chr == "CN",], "black", highlightSexChrs, pch=15, cex=cex.points+0.5, highlightCnPoints=highlightCnPoints)
addPoints(data[data$chr == "CN",], "yellow", highlightSexChrs, pch=15, cex=cex.points, highlightCnPoints=highlightCnPoints)
}
##TODO - add a legend for point types - highlight vs sex chrs vs autosomes
## add highlighted points selected for by user
if(!(is.null(positionsToHighlight))){
addpts = merge(data, positionsToHighlight, by.x=c("chr","st"), by.y = c("chr","st"))
xs <- c()
ys <- c()
labels <- c()
if(length(addpts[,1]) > 0){
if(highlightsHaveNames){
if(labelOnPlot) {
xs <- addpts$vaf[addpts$name != ""]
ys <- addpts$depth[addpts$name != ""]
labels <- addpts$name[addpts$name != ""]
num.labels <- length(labels)
suppressPackageStartupMessages(library(TeachingDemos))
# By using the min argument, ensure that the annotations
# do not overlap the corresponding symbol. NB: we anticipate
# this code will only be active for the interior points
text(x=addpts$vaf,y=addpts$depth,label=rep("*", length(addpts$vaf)))
if(num.labels > 0){
xs <- spread.labs(xs, mindiff=4, min=xs)
ys <- spread.labs(ys+50, mindiff=4, min=ys)
text(x=xs,y=ys,label=labels)
}
} else {
for(i in 1:length(addpts$vaf)){
if(addpts$name[i] != "") {
# Only label the gene with a number if it has a number
text(addpts$vaf[i],addpts$depth[i],labels=i,cex=0.5)
} else {
# Otherwise, just highlight it with a star
text(addpts$vaf[i],addpts$depth[i],labels="*")
}
}
}
} else {
addPoints(addpts, col="#555555FF", highlightSexChrs, cex=cex.points, highlightCnPoints=highlightCnPoints);
}
}
}
}
## define the axis
axis(side=2,las=1,tck=0,lwd=0,cex.axis=1.2/(scale*2)*textScale, hadj=0.5/scale)
for (i in 2:length(axTicks(2)-1)) {
lines(c(-1,101),c(axTicks(2)[i],axTicks(2)[i]),col="#00000022");
}
## plot the background color
rect(-1, 5, 101, axTicks(2)[length(axTicks(2))]*1.05, col = "#00000011",border=NA);
## add cn circle
cnpos = axTicks(2)[length(axTicks(2))-1]
points(x=c(95), y=cnpos, type="p", pch=19, cex=3/scale, col=circlecolor);
text(c(95),y=cnpos, labels=c(cn), cex=1/scale, col="#FFFFFFFF")
## y axis label
mtext("Tumor Coverage", side=2, cex=0.6*(textScale*0.8), padj=-4.2*(axisPosScale));
}
##--------------------------------------------------------------------
## add a legend for highlighted points with names
addHighlightLegend <- 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() }
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);
# mtext("Genes",side=2,cex=0.5,padj=-4.2);
ypos=rev(seq(0,900,(900/13)))[1:13]
# Don't put empty names on legend
non.trivial.indices <- (1:length(addpts$name))[addpts$name != ""]
ncol=ceiling(length(non.trivial.names)/13)
xpos=0;
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(non.trivial.indices[nxt],". ",names[i],sep=""), cex=0.6, pos=4)
text(xpos, ypos[i], paste(nxt,". ",names[i],sep=""), cex=0.6, pos=4)
nxt <- nxt+1
}
xpos=xpos+250;
offset=offset+13;
}
}
##---------------------------------------------------------------------
## create the two dimensional plot with scatter annotated with
## clustering results and 1D plots along margins, this time using
## ggplot2
sc.plot2dWithMargins <- function(sco, outputFile, positionsToHighlight=NULL, highlightsHaveNames=FALSE, overlayErrorBars=FALSE) {
densityData = sco@densities
vafs.merged = sco@vafs.merged
vafs.1d.merged = sco@vafs.1d
sampleNames = sco@sampleNames
dimensions = sco@dimensions
clust = sco@clust
marginalClust = sco@marginalClust
suppressPackageStartupMessages(library(grid))
suppressPackageStartupMessages(library(ggplot2))
plots.1d.list <- list()
res.1d.max.densities <- list()
xmin <- -5
xmax <- 105
# NB: ymin and ymax below on the right-hand side of the = are
# referring to verbatim names in the data argument. But R CMD check
# doesn't understand this, so define these variables. This is only
# a hack to quiet R CMD check and does not affect the ggplot call.
ymin <- NULL
ymax <- NULL
tmp.file <- tempfile("outputFile.tmp")
pdf(file=tmp.file, width=7.2, height=6, bg="white")
# Create (and store) all possible 1D plots
for(d in 1:dimensions){
# Overlay the histogram of the data on the model fit.
#ylab <- "\nDensity\n"
ylab <- "Density (a.u.)"
title <- ""
# Set max posterior density to max of splinefun
limits <- data.frame(x=c(min(marginalClust[[d]]$fit.x), max(marginalClust[[d]]$fit.x)))
f <- splinefun(marginalClust[[d]]$fit.x, marginalClust[[d]]$fit.y[1,,drop=FALSE])
max.posterior.density <- max(unlist(lapply(seq(from=limits$x[1], to=limits$x[2], by=10^-3), f)))
# xlab <- paste("\n", sampleNames[d], "\n", sep="")
xlab <- paste(sampleNames[d], " VAF", sep="")
numClusters = 0
if(!is.null(vafs.1d.merged[[d]]$cluster)) {
numClusters = max(vafs.1d.merged[[d]]$cluster, na.rm=T)
}
vafs = getOneSampleVafs(vafs.1d.merged[[d]], sampleNames[d], numClusters)
# Only show copy number = 2
v = vafs[which(vafs$cleancn==2 & vafs$adequateDepth==1),];
frequencies <- data.frame(x=v$vaf, row.names=NULL, stringsAsFactors=NULL)
g <- ggplot(data = frequencies, aes(x)) + ggtitle(title) + xlab(xlab) + ylab(ylab)
# g <- g + theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank())
g <- g + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank())
g <- g + theme(plot.margin = unit(c(0,0,0,0), "cm"))
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)))
# This is just to appease R CMD check--it does not effect the
# geom_histogram call below.
..ncount..=NULL;
g <- g + geom_histogram(data = frequencies, mapping=aes(x, y=..ncount..*100), fill="white", colour=gray(0.4), breaks=breaks)
# Need to "print" the graph in order to see its maximum y value
# NB: this is redundant at this point, given that I scale the
# histogram to have height 100 above. i.e., max.density will be 100.
tmp <- print(g)
max.density <- max(tmp[["data"]][[1]]$ymax)
res.1d.max.densities[[d]] <- max.density
hline <- data.frame(x = c(0,100), y=c(-5,-5))
g <- g + geom_line(data = hline, aes(x,y))
#vline <- data.frame(y = c(0,max.density * 1.1), x=c(-5,-5))
vline <- data.frame(y = c(0,100), x=c(-5,-5))
g <- g + geom_line(data = vline, aes(x,y))
scale <- max.density / max.posterior.density
# If we actually performed clustering, show it.
if(numClusters > 0) {
f <- splinefun(marginalClust[[d]]$fit.x, scale*marginalClust[[d]]$fit.y[1,,drop=FALSE])
g <- g + stat_function(data = limits, fun=f, mapping=aes(x))
}
g <- g + coord_cartesian(ylim=c(xmin, max.density*1.1), xlim=c(xmin,xmax))
plots.1d.list[[d]] <- g
}
devoff = dev.off()
unlink(tmp.file)
vplayout <- function(x, y) viewport(layout.pos.row = x, layout.pos.col = y, clip="off")
pdf(file=outputFile, width=7.2, height=6, bg="white")
##create a 2d plot for each pairwise combination of samples
for(d1 in 1:(dimensions-1)){
for(d2 in d1:dimensions){
if(d1==d2){
next
}
xlab <- sampleNames[d1]
xlab <- gsub("\\.", " ", xlab)
#xlab <- paste("\n", xlab, "\n", sep="")
xlab <- paste(xlab, " VAF", sep="")
ylab <- sampleNames[d2]
ylab <- gsub("\\.", " ", ylab)
#ylab <- paste("\n", ylab, "\n", sep="")
ylab <- paste(ylab, " VAF", sep="")
numClusters = 0
if(!is.null(vafs.merged$cluster)) {
numClusters = max(vafs.merged$cluster, na.rm=T)
}
vafs1 = getOneSampleVafs(vafs.merged, sampleNames[d1], numClusters);
vafs2 = getOneSampleVafs(vafs.merged, sampleNames[d2], numClusters);
##get only cn2 points
vafs1 = vafs1[which(vafs1$cleancn==2 & vafs1$adequateDepth==1),]
vafs2 = vafs2[which(vafs2$cleancn==2 & vafs2$adequateDepth==1),]
if(numClusters > 0) {
v = merge(vafs1,vafs2,by.x=c("chr","st","cluster"), by.y=c("chr","st","cluster"),suffixes=c(".1",".2"))
# Remove any outliers--these will have cluster assignment 0
v.outlier <- v[v$cluster == 0,]
v <- v[v$cluster != 0,]
} else {
v = merge(vafs1,vafs2,by.x=c("chr","st"), by.y=c("chr","st"),suffixes=c(".1",".2"))
}
v.no.highlight <- v
if(!(is.null(positionsToHighlight))) {
names(positionsToHighlight)=c("chr","st","name");
chr.start.v <- cbind(v[,"chr"], v[,"st"])
chr.start.highlight <- cbind(positionsToHighlight[,1], positionsToHighlight[,2])
v.no.highlight <- v[!(apply(chr.start.v, 1, paste, collapse="$$") %in% apply(chr.start.highlight, 1, paste, collapse="$$")),]
}
title <- ""
# Plot the points that we will not highlight
frequencies.no.highlight <- data.frame(x=v.no.highlight$vaf.1, y=v.no.highlight$vaf.2, row.names=NULL, stringsAsFactors=NULL)
colvec = c()
if(numClusters > 0) {
clusters <- v.no.highlight$cluster
cols=getClusterColors(numClusters)
colvec = cols[clusters]
clusters <- unlist(lapply(clusters, as.character))
frequencies.no.highlight <- data.frame(x=v.no.highlight$vaf.1, y=v.no.highlight$vaf.2, shape=clusters, colour=clusters, row.names=NULL, stringsAsFactors=FALSE)
# NB: shape and colour below on the right-hand side of the = are
# referring to verbatim names in the data argument. But R CMD check
# doesn't understand this, so define these variables. This is only
# a hack to quiet R CMD check and does not affect the ggplot call.
shape <- NULL
colour <- NULL
g <- ggplot(data = frequencies.no.highlight, aes(x=x, y=y)) + ggtitle(title) + xlab(xlab) + ylab(ylab) + geom_point(data = frequencies.no.highlight, aes(x=x, y=y, shape=shape, colour=colour))
# Plot a legend
g <- g + scale_colour_manual(name = "Clusters", labels = 1:numClusters, values = cols[1:numClusters])
g <- g + scale_shape_manual(name = "Clusters", labels = 1:numClusters, values = 1:numClusters)
} else {
g <- ggplot(data = frequencies.no.highlight, aes(x=x, y=y)) + ggtitle(title) + xlab(xlab) + ylab(ylab) + geom_point(data = frequencies.no.highlight, aes(x=x, y=y))
}
if(overlayErrorBars == TRUE) {
err.bars.1 <- compute.binomial.error.bars(v.no.highlight$var.1, v.no.highlight$depth.1) * 100
err.bars.2 <- compute.binomial.error.bars(v.no.highlight$var.2, v.no.highlight$depth.2) * 100
err.df.x <- data.frame(x=v.no.highlight$vaf.1, y=v.no.highlight$vaf.2, xmin=err.bars.1$lb, xmax=err.bars.1$ub)
err.df.y <- data.frame(x=v.no.highlight$vaf.1, y=v.no.highlight$vaf.2, ymin=err.bars.2$lb, ymax=err.bars.2$ub)
if(numClusters > 0) {
g <- g + geom_errorbarh(data = err.df.x, aes(x=x, y=y, xmin=xmin, xmax=xmax), colour=colvec)
g <- g + geom_errorbar(data = err.df.y, aes(x=x, y=y, ymin=ymin, ymax=ymax), colour=colvec)
} else {
g <- g + geom_errorbarh(data = err.df.x, aes(x=x, y=y, xmin=xmin, xmax=xmax))
g <- g + geom_errorbar(data = err.df.y, aes(x=x, y=y, ymin=ymin, ymax=ymax))
}
}
# Now overlay any points that we will highlight
if(!(is.null(positionsToHighlight))) {
# Merge the data and the positions to highlight by chr (col 1)
# and start (col 2)
addpts = merge(v, positionsToHighlight, by.x=c("chr","st"), by.y = c("chr","st"))
if(dim(addpts)[1] > 0) {
frequencies.highlight <- data.frame(x=addpts$vaf.1, y=addpts$vaf.2, row.names=NULL, stringsAsFactors=NULL)
g <- g + geom_point(data = frequencies.highlight, aes(x=x, y=y), shape="*", size=10, colour="black")
if(overlayErrorBars == TRUE) {
err.bars.1 <- compute.binomial.error.bars(addpts$var.1, addpts$depth.1) * 100
err.bars.2 <- compute.binomial.error.bars(addpts$var.2, addpts$depth.2) * 100
err.df.x <- data.frame(x=addpts$vaf.1, y=addpts$vaf.2, xmin=err.bars.1$lb, xmax=err.bars.1$ub)
g <- g + geom_errorbarh(data = err.df.x, aes(x=x, y=y, xmin=xmin, xmax=xmax), colour="black")
err.df.y <- data.frame(x=addpts$vaf.1, y=addpts$vaf.2, ymin=err.bars.2$lb, ymax=err.bars.2$ub)
g <- g + geom_errorbar(data = err.df.y, aes(x=x, y=y, ymin=ymin, ymax=ymax), colour="black")
}
}
}
# g <- g + theme_bw() + theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank())
g <- g + theme_bw() + theme(panel.border = element_blank())
g <- g + theme(plot.margin = unit(c(0,0,0,0), "cm"))
# Put the legend in the top right
g <- g + theme(legend.position = c(1,1), legend.justification=c(1,1))
hline <- data.frame(x = c(0,100), y=c(-5,-5))
g <- g + geom_line(data = hline, aes(x,y))
vline <- data.frame(y = c(0,100), x=c(-5,-5))
g <- g + geom_line(data = vline, aes(x,y))
g <- g + coord_cartesian(xlim=c(xmin, xmax), ylim=c(xmin, xmax))
plot.2d <- g
# Code to override clipping
#plot.2d <- ggplot_gtable(ggplot_build(g))
#plot.2d$layout$clip[plot.2d$layout$name=="panel"] <- "off"
#grid.draw(plot.2d)
plot.1d.1 <- plots.1d.list[[d1]]
plot.1d.2 <- plots.1d.list[[d2]]
grid.newpage()
pushViewport(viewport(layout = grid.layout(2, 2), clip="off"))
text.size <- 10
vp <- vplayout(1,1)
vp11 <- vp
plot.2d <- plot.2d + theme(text = element_text(size = text.size))
print(plot.2d, vp=vp)
if(!(is.null(positionsToHighlight))) {
# Merge the data and the positions to highlight by chr (col 1)
# and start (col 2)
addpts = merge(v, positionsToHighlight, by.x=c("chr","st"), by.y = c("chr","st"))
frequencies.highlight <- data.frame(x=addpts$vaf.1, y=addpts$vaf.2, row.names=NULL, stringsAsFactors=NULL)
# Add the labels
if(dim(addpts)[1] > 0){
if(highlightsHaveNames){
# This code adapted from stackoverflow.com/questions/10536396/using-grconvertx-grconverty-in-ggplot2
# Create a new viewport with clipping disabled so we can
# put text outside the plot
depth <- downViewport('panel.3-4-3-4')
pushViewport(dataViewport(xData=c(0,100), yData=c(0,100), clip='off'))
xs <- list()
ys <- list()
labels <- list()
nxt <- 1
for(i in 1:dim(addpts)[1]) {
if(addpts$vaf.1[i] < 1) {
x <- addpts$vaf.1[i] - 12
y <- addpts$vaf.2[i]
} else if(addpts$vaf.2[i] < 1) {
x <- addpts$vaf.1[i]
y <- addpts$vaf.2[i] - 10
} else {
x <- addpts$vaf.1[i] + 5
y <- addpts$vaf.2[i] + 5
}
label <- as.character(addpts$name[i])
#df <- data.frame(x=x, y=y)
if(label == "") {
# No label to place
} else {
# grid.text(x=x,y=y,label=label,default.units="native", gp=gpar(fontsize=8))
xs[[nxt]] <- x
ys[[nxt]] <- y
labels[[nxt]] <- label
nxt <- nxt + 1
}
}
xs <- unlist(xs)
ys <- unlist(ys)
labels <- unlist(labels)
num.labels <- length(labels)
suppressPackageStartupMessages(library(TeachingDemos))
# By using the min argument, ensure that the annotations
# do not overlap the corresponding symbol. NB: we anticipate
# this code will only be active for the interior points
if(num.labels > 0){
xs <- spread.labs(xs, mindiff=4, min=xs)
ys <- spread.labs(ys, mindiff=4, min=ys)
for(i in 1:num.labels) {
grid.text(x=xs[i],y=ys[i],label=labels[i],default.units="native", gp=gpar(fontsize=8))
}
}
# Move up depth+1 in the tree. NB: +1 because we pushed a
# viewport on to the tree.
upViewport(depth+1)
}
}
}
vp <- vplayout(1,2)
plot.1d.2 <- plot.1d.2 + theme(text = element_text(size = text.size))
plot.1d.2 <- plot.1d.2 + coord_flip(ylim = c(xmin, res.1d.max.densities[[d2]]*1.1), xlim = c(xmin, xmax))
print(plot.1d.2, vp=vp)
vp <- vplayout(2,1)
plot.1d.1 <- plot.1d.1 + theme(text = element_text(size = text.size))
print(plot.1d.1, vp=vp)
}
}
devoff = dev.off()
}
## Compute the lower and upper bound for a "1-std dev" binomial confidence
## interval for a given number of successes and total number trials
compute.binomial.error.bars <- function(successes, total.trials){
suppressPackageStartupMessages(library(MKmisc))
suppressPackageStartupMessages(library(NORMT3))
# Return a "1 std dev" confidence interval
width <- as.double(erf(1/sqrt(2)))
lb <- mapply(function(a,b) binomCI(a, b, conf.level=width, method="jeffreys")$CI[1], successes, total.trials)
ub <- mapply(function(a,b) binomCI(a, b, conf.level=width, method="jeffreys")$CI[2], successes, total.trials)
res <- data.frame(lb=lb, ub=ub)
return(res)
}
##---------------------------------------------------------------------
## create the two dimensional plot with scatter annotated with
## clustering results and 1D plots along margins
##---------------------------------------------------------------------------------
## Create two dimensional plot with scatter annotated with clustering result
##
sc.plot2d <- function(sco, outputFile=NULL, positionsToHighlight=NULL, highlightsHaveNames=FALSE, overlayClusters=TRUE, overlayErrorBars=FALSE, ellipse.metadata = list(), singlePage=FALSE, scale=1, xlim=100, ylim=100, plot.title=NULL, samplesToPlot=NULL, clusterLegend=TRUE, flipSamples=FALSE, xlab=NULL, ylab=NULL, colors=NULL, plotWidth=7.2, plotHeight=6){
vafs.merged = sco@vafs.merged
sampleNames = sco@sampleNames
dimensions = sco@dimensions
if(singlePage){
nplots = ncol(combn(c(1:dimensions),2))
nrow=round(sqrt(nplots))
ncol=ceiling(sqrt(nplots))
if(!is.null(outputFile)){
pdf(outputFile, width=plotWidth*ncol, height=plotHeight*nrow, bg="white")
}
par(mfrow=c(nrow,ncol), mar=c(5.1, 5.1, 4.1, 2.1))
} else {
if(!is.null(outputFile)){
pdf(outputFile, width=plotWidth, height=plotHeight, bg="white")
}
}
numClusters = 0
maxCluster = 0
non.empty.clusters <- c()
# NB: clusters no longer need be numbered contiguously
if(!is.null(vafs.merged$cluster)) {
# Careful! This doesn't work--0 is an outlier cluster.
#numClusters = length(unique(vafs.merged$cluster[!is.na(vafs.merged$cluster)]))
maxCluster = max(vafs.merged$cluster, na.rm=T)
# Note that by starting at 1, we exclude the outlier "cluster" 0.
for(i in 1:maxCluster){
if(!any(vafs.merged$cluster == i)) { next }
non.empty.clusters <- c(non.empty.clusters, i)
}
}
numClusters = length(non.empty.clusters)
suppressPackageStartupMessages(library(plotrix)) # For plotCI among others.
createPlot <- function(d1, d2){
vafs1 = getOneSampleVafs(vafs.merged, sampleNames[d1], numClusters);
vafs2 = getOneSampleVafs(vafs.merged, sampleNames[d2], numClusters);
##get only cn2 points with adequate coverage
vafs1 = vafs1[which(vafs1$cleancn==2 & vafs1$adequateDepth==1),]
vafs2 = vafs2[which(vafs2$cleancn==2 & vafs2$adequateDepth==1),]
if(!is.null(vafs.merged$cluster)) {
v = merge(vafs1,vafs2,by.x=c("chr","st","cluster"), by.y=c("chr","st","cluster"),suffixes=c(".1",".2"))
## Remove any outliers--these will have cluster assignment 0
v.outlier <- v[v$cluster == 0,]
v <- v[v$cluster != 0,]
} else {
v = merge(vafs1,vafs2,by.x=c("chr","st"), by.y=c("chr","st"),suffixes=c(".1",".2"))
}
cols=colors
if(is.null(colors)){
cols = getClusterColors(maxCluster)
}
##sample name
title=paste(sampleNames[d1],"vs",sampleNames[d2])
if(!is.null(plot.title)){
title=plot.title
}
##create the plot
xlabel = paste(sampleNames[d1],"VAF ")
ylabel = paste(sampleNames[d2],"VAF")
if(!(is.null(xlab))){
xlabel = xlab
}
if(!(is.null(ylab))){
ylabel = ylab
}
plot(-100, -100, xlim=c(0,xlim*1.2), ylim=c(0,ylim), main=title,
xlab=xlabel, ylab=ylabel,
bty="n", xaxt="n", yaxt="n", cex.lab=scale, cex.main=scale, cex.axis=scale)
xGridIncrement=20
yGridIncrement=20
if(xlim < 80){
xGridIncrement= 10
}
if(xlim < 50){
xGridIncrement=5
}
if(ylim < 80){
yGridIncrement=10
}
if(ylim < 50){
yGridIncrement=5
}
## vertical grid
abline(v=seq(0,xlim,xGridIncrement),col="grey50", lty=3, lwd=scale)
axis(side=1,at=seq(0,xlim,xGridIncrement),labels=seq(0,xlim,xGridIncrement), cex.axis=scale)
## horizontal grid
numYGridLines=floor(ylim/yGridIncrement)
segments(rep(-10,numYGridLines),seq(0,ylim,yGridIncrement),rep(xlim*1.05,numYGridLines),seq(0,ylim,yGridIncrement), lty=3, lwd=scale, col="grey50")
axis(side=2,at=seq(0,ylim,yGridIncrement),labels=seq(0,ylim,yGridIncrement), cex.axis=scale)
## If we will be highlighting some points, exclude them from
## the general list of points to plot and plot them instead with
## a different symbol/color (a black *)
v.no.highlight <- v
if(!(is.null(positionsToHighlight))) {
names(positionsToHighlight)=c("chr","st","name");
chr.start.v <- cbind(v[,"chr"], v[,"st"])
chr.start.highlight <- cbind(positionsToHighlight[,1], positionsToHighlight[,2])
v.no.highlight <- v[!(apply(chr.start.v, 1, paste, collapse="$$") %in% apply(chr.start.highlight, 1, paste, collapse="$$")),]
}
if(overlayErrorBars == TRUE) {
err.bars.1 <- compute.binomial.error.bars(v.no.highlight$var.1, v.no.highlight$depth.1) * 100
err.bars.2 <- compute.binomial.error.bars(v.no.highlight$var.2, v.no.highlight$depth.2) * 100
}
if(!is.null(vafs.merged$cluster)) {
## handle cluster 0 (outliers)
if(length(v.outlier[,1]) > 0){
points(v.outlier$vaf.1, v.outlier$vaf.2, col=rgb(0,0,0,0.5), pch=".", cex=3*scale)
}
for(i in 1:numClusters){
indices <- v.no.highlight$cluster==non.empty.clusters[i]
if(overlayClusters){
if(dim(v.no.highlight[indices,])[1] > 0) {
if(overlayErrorBars == TRUE) {
plotCI(v.no.highlight[indices,]$vaf.1, v.no.highlight[indices,]$vaf.2, col=cols[non.empty.clusters[i]], pch=i, li=err.bars.1$lb[indices], ui=err.bars.1$ub[indices], add=TRUE, err="x")
plotCI(v.no.highlight[indices,]$vaf.1, v.no.highlight[indices,]$vaf.2, col=cols[non.empty.clusters[i]], pch=i, li=err.bars.2$lb[indices], ui=err.bars.2$ub[indices], add=TRUE, err="y")
} else {
points(v.no.highlight[indices,]$vaf.1, v.no.highlight[indices,]$vaf.2, col=cols[non.empty.clusters[i]], pch=i, cex=scale)
}
}
} else { ##no overlay of clusters
if(dim(v.no.highlight[indices,])[1] > 0) {
if(overlayErrorBars == TRUE) {
plotCI(v.no.highlight[indices,]$vaf.1, v.no.highlight[indices,]$vaf.2, pch=18, li=err.bars.1$lb[indices], ui=err.bars.1$ub[indices], add=TRUE, err="x")
plotCI(v.no.highlight[indices,]$vaf.1, v.no.highlight[indices,]$vaf.2, pch=18, li=err.bars.2$lb[indices], ui=err.bars.2$ub[indices], add=TRUE, err="y")
} else {
points(v.no.highlight[indices,]$vaf.1, v.no.highlight[indices,]$vaf.2, pch=18, cex=scale,col=getClusterColors(1)[1])
}
}
}
}
} else { #no clusters plotted
if(dim(v.no.highlight)[1] > 0) {
if(overlayErrorBars == TRUE) {
plotCI(v.no.highlight$vaf.1, v.no.highlight$vaf.2, pch=18, li=err.bars.1$lb, ui=err.bars.1$ub, add=TRUE, err="x", cex=scale, col=getClusterColors(1)[1])
plotCI(v.no.highlight$vaf.1, v.no.highlight$vaf.2, pch=18, li=err.bars.2$lb, ui=err.bars.2$ub, add=TRUE, err="y", cex=scale, col=getClusterColors(1)[1])
} else {
points(v.no.highlight$vaf.1, v.no.highlight$vaf.2, pch=18, cex=scale, col=getClusterColors(1)[1])
}
}
}
## Now plot the highlighted points so they are overlaid
if(!(is.null(positionsToHighlight))) {
## Merge the data and the positions to highlight by chr (col 1)
## and start (col 2)
addpts = merge(v, positionsToHighlight, by.x=c("chr","st"), by.y = c("chr","st"))
## Plot the highlighted items. NB: we never overlay the
## cluster on them, but expect this will be obvious from context
if(dim(addpts)[1] > 0) {
if(overlayErrorBars == TRUE) {
err.bars.1 <- compute.binomial.error.bars(addpts$var.1, addpts$depth.1) * 100
err.bars.2 <- compute.binomial.error.bars(addpts$var.2, addpts$depth.2) * 100
plotCI(addpts$vaf.1, addpts$vaf.2, pch="*", col="black", cex=2, li=err.bars.1$lb, ui=err.bars.1$ub, add=TRUE, err="x")
plotCI(addpts$vaf.1, addpts$vaf.2, pch="*", col="black", cex=2, li=err.bars.2$lb, ui=err.bars.2$ub, add=TRUE, err="y")
} else {
points(addpts$vaf.1, addpts$vaf.2, pch="*", col="black", cex=2*scale)
}
}
}
if(!is.null(vafs.merged$cluster)) {
for(i in 1:numClusters){
if((!is.null(ellipse.metadata$SEMs.lb)) & (!is.null(ellipse.metadata$SEMs.ub))) {
xc <- ellipse.metadata$SEMs.lb[i,d1] + ((ellipse.metadata$SEMs.ub[i,d1] - ellipse.metadata$SEMs.lb[i,d1])/2)
yc <- ellipse.metadata$SEMs.lb[i,d2] + ((ellipse.metadata$SEMs.ub[i,d2] - ellipse.metadata$SEMs.lb[i,d2])/2)
## ell <- my.ellipse(hlaxa = ((ellipse.metadata$SEMs.ub[i,d1] - ellipse.metadata$SEMs.lb[i,d1])/2), hlaxb = ((ellipse.metadata$SEMs.ub[i,d2] - ellipse.metadata$SEMs.lb[i,d2])/2), xc = xc, yc = yc)
draw.ellipse(xc, yc, a = ((ellipse.metadata$SEMs.ub[i,d1] - ellipse.metadata$SEMs.lb[i,d1])/2), b = ((ellipse.metadata$SEMs.ub[i,d2] - ellipse.metadata$SEMs.lb[i,d2])/2))
}
if((!is.null(ellipse.metadata$std.dev.lb)) & (!is.null(ellipse.metadata$std.dev.ub))) {
xc <- ellipse.metadata$std.dev.lb[i,d1] + ((ellipse.metadata$std.dev.ub[i,d1] - ellipse.metadata$std.dev.lb[i,d1])/2)
yc <- ellipse.metadata$std.dev.lb[i,d2] + ((ellipse.metadata$std.dev.ub[i,d2] - ellipse.metadata$std.dev.lb[i,d2])/2)
## Plot std dev as dashed line.
draw.ellipse(xc, yc, a = ((ellipse.metadata$std.dev.ub[i,d1] - ellipse.metadata$std.dev.lb[i,d1])/2), b = ((ellipse.metadata$std.dev.ub[i,d2] - ellipse.metadata$std.dev.lb[i,d2])/2), lty=2)
}
}
}
##plot(-100, -100, xlim=c(0,100), ylim=c(0,100),main="Clusters")
if(!is.null(vafs.merged$cluster)) {
## Only include in the legend any clusters that are not empty.
if(clusterLegend==TRUE){
legend("topright", legend=non.empty.clusters, col=cols[non.empty.clusters], title="Clusters", pch=non.empty.clusters, cex=scale, bg="white")
}
##legend("topright", legend=1:numClusters, col=cols[1:numClusters], title="Clusters", pch=1:numClusters)
}
## Add annotation for gene names, if requested
if(!(is.null(positionsToHighlight))){
addpts = merge(v, positionsToHighlight, by.x=c("chr","st"), by.y = c("chr","st"))
## write.table(file="genes.txt", unique(addpts$gene_name), row.names=FALSE, col.names=FALSE, quote=FALSE)
if(dim(addpts)[1] > 0){
if(highlightsHaveNames){
xs <- list()
ys <- list()
labels <- list()
nxt <- 1
for(i in 1:dim(addpts)[1]) {
par(xpd=NA)
cex <- 1
if(addpts$vaf.1[i] < 1) {
##text(addpts$vaf.1[i] - 8,addpts$vaf.2[i],labels=addpts$gene_name[i],cex=cex)
x <- addpts$vaf.1[i] - 8
y <- addpts$vaf.2[i]
} else if(addpts$vaf.2[i] < 1) {
##text(addpts$vaf.1[i],addpts$vaf.2[i] - 5,labels=addpts$gene_name[i],cex=cex)
x <- addpts$vaf.1[i]
y <- addpts$vaf.2[i] - 5
} else {
##text(addpts$vaf.1[i] + 4,addpts$vaf.2[i] + 4,labels=addpts$gene_name[i],cex=cex)
x <- addpts$vaf.1[i] + 4
y <- addpts$vaf.2[i] + 4
}
label <- as.character(addpts$name[i])
##df <- data.frame(x=x, y=y)
##grid.text(x=x,y=y,label=label,default.units="native", gp=gpar(fontsize=8))
if(label != "") {
xs[[nxt]] <- x
ys[[nxt]] <- y
labels[[nxt]] <- label
nxt <- nxt + 1
}
}
xs <- unlist(xs)
ys <- unlist(ys)
labels <- unlist(labels)
num.labels <- length(labels)
suppressPackageStartupMessages(library(TeachingDemos))
## By using the min argument, ensure that the annotations
## do not overlap the corresponding symbol. NB: we anticipate
## this code will only be active for the interior points
if(num.labels > 0){
xs <- spread.labs(xs, mindiff=4, min=xs)
ys <- spread.labs(ys, mindiff=4, min=ys)
for(i in 1:num.labels) {
##grid.text(x=xs[i],y=ys[i],label=labels[i],default.units="native", gp=gpar(fontsize=8))
text(x=xs[i], y=ys[i], label=labels[i], cex=cex)
}
}
}
}
} # End add gene annotations
}
##create 2d plots
count=0;
for(d1 in 1:(dimensions-1)){
for(d2 in d1:dimensions){
if(d1==d2){
next
}
if(is.null(samplesToPlot)){ ##use each pairwise combination of samples
if(flipSamples == TRUE){
createPlot(d2,d1)
} else {
createPlot(d1,d2)
}
count = count + 1;
} else { ##just plot the specified samples
if( (sampleNames[d1] %in% samplesToPlot) & (sampleNames[d2] %in% samplesToPlot)){
if(flipSamples == TRUE){
createPlot(d2,d1)
} else {
createPlot(d1,d2)
}
count = count + 1;
}
}
}
}
if(count < 1){
print ("WARNING: no pair of samples matched the samplesToPlot argument. Nothing plotted");
}
if(!is.null(outputFile)){
devoff = dev.off()
}
}
##-------------------------------------------------------------------------------------
## plot three samples in 3d, optionally create a GIF
##
sc.plot3d <- function(sco, samplesToPlot, size=700, outputFile=NULL){
suppressPackageStartupMessages(library(rgl))
##set the size of the window
open3d(windowRect = c(0,0, size, size) )
if(length(samplesToPlot) != 3){
print("ERROR: must provide exactly 3 sample names for 3d plotting")
return(0)
}
a = sco@vafs.merged[,c(paste(samplesToPlot,".vaf",sep=""),"cluster")]
a = a[!is.na(a$cluster),]
a = a[!(a$cluster==0),]
numClusters=max(a$cluster)
cols=getClusterColors(numClusters)
colvec = cols[a$cluster]
plot3d(a[,1], a[,2], a[,3], xlim=c(0,100), ylim=c(0,100),zlim=c(0,100), axes=FALSE,
xlab=samplesToPlot[1], ylab=samplesToPlot[2], zlab=samplesToPlot[3],
type="s", col=colvec)
##add a box
axes3d( edges=c("x--", "y--", "z"),labels=FALSE)
for(i in c("+","-")){
for(j in c("+","-")){
axes3d( edges=paste("x",i,j,sep=""), tick=FALSE, labels=FALSE)
axes3d( edges=paste("y",i,j,sep=""), tick=FALSE, labels=FALSE)
axes3d( edges=paste("z",i,j,sep=""), tick=FALSE, labels=FALSE)
}
}
if(is.null(outputFile)){
play3d(spin3d(axis=c(0,0,1), rpm=10), duration=6)
} else {
##remove trailing .gif, since movie3d adds it
outputFile = sub(".gif$","",outputFile)
movie3d(spin3d(axis=c(0,0,1), rpm=5), duration=12, dir=getwd(), movie=outputFile)
}
}
##-------------------------------------------------------------------------------------
## Get a list of colors to use for the clusters (after 20 they start repeating)
##
getClusterColors <- function(numClusters){
suppressPackageStartupMessages(library(RColorBrewer))
cols=brewer.pal(8,"Dark2")
if(numClusters>8){
cols = c(cols,brewer.pal(9,"Set1"))
}
if(numClusters > 20){ #if this happens, something is probably weird with your data, but we'll support it anyway
cols = rep(cols,ceiling(length(cols)/20))
}
cols = cols[1:numClusters]
##add transparency to colors
for(i in 1:length(cols)){
z = col2rgb(cols[i])
cols[i] = rgb(z[1], z[2], z[3], 150, maxColorValue=255)
}
return(cols)
}
##-------------------------------------------------------------------------------------
## Extract a single sample's data from the merged data frame
##
getOneSampleVafs <- function(vafs.merged, name, numClusters){
common = c("chr","st","adequateDepth")
a = vafs.merged[,common]
prefix = paste("^",name,".",sep="")
cols=grepl(prefix, names(vafs.merged))
header=sub(prefix, "", names(vafs.merged)[cols])
vafs = cbind(a,vafs.merged[,cols])
names(vafs) = c(common,header)
if(numClusters > 0){
vafs = cbind(vafs, vafs.merged$cluster)
names(vafs)[length(vafs)] = "cluster"
}
return(vafs)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.