## Note: A simple R script to produce high quality figure of mutational landscape
##+for paper, developed and maintained by Li Xiangchun in BGI, all rights reserved.
## It's highly recommended NOT to include '-' in sample names.
plot.land <- function(upper.panel, middle.panel, left.panel, right.panel, legend.panel, subtype.panel=NULL, subtype.2.panel=NULL,
output.fig.pdffile=NULL,
left.panel.xlim=NULL,
family.font="sans",
syn.nonsyn.legend.color=NULL,
syn.nonsyn.legend.ncol=NULL,
syn.nonsyn.legend.loc="top",
image.color=NULL,
subtype.labels=NULL,
subtype.label.cex=1,
subtype.panel.color=brewer.pal(7, "Dark2"),
subtype.panel.identifier="",
subtype.panel.identifier.cex=0.6,
subtype.legend.text.width=NULL,
subtype.2.labels=NULL,
subtype.2.label.cex=1,
subtype.2.panel.color=brewer.pal(7, "Set2"),
subtype.2.panel.identifier="",
subtype.2.panel.identifier.cex=0.6,
subtype.2.legend.text.width=NULL,
right.panel.color=brewer.pal(7, "Purples")[5],
sort.by.qvalue=FALSE,
group.by.subtype=FALSE,
left.panel.cex.axis=1.5,
right.panel.cex.axis=1.3,
upper.panel.mar=c(0, 1.3, 0.5, 1.16),
legend.panel.mar=c(0, 1, 2, 1),
legend.panel.ncol=NULL,
legend.panel.y=-0.05,
left.panel.mar=c(5.5, 1.5, 0.2, 2.9),
middle.panel.mar=c(5.5, 1.3, 0.1, 1.16),
right.panel.mar=c(5.5, 3, 0.2, 1.5),
left.panel.xlab="Individuals\nwith mutations",
left.panel.xlab.cex=0.8,
right.panel.xlab="-log10\n(q-value)",
right.panel.xlab.cex=0.8,
right.panel.qscore.cutoff=1,
right.panel.qscore.lwd=1,
right.panel.qscore.lty='dashed',
right.panel.qscore.col='red')
{
upper.panel.samples <- as.character(upper.panel[,1])
middle.panel.samples <- as.character(colnames(middle.panel))
left.panel.genes <- as.character(left.panel[,1])
right.panel.genes <- as.character(right.panel[,1])
if (length(upper.panel.samples) != length(middle.panel.samples)) {
stop("The number of samples in upper panel is NOT equal to the number of samples in middle panel!")
}
if (!all(sort(upper.panel.samples) == sort(middle.panel.samples))) {
stop("Sample names in upper panel are not all of the same as those in middle panel!")
}
# Some constant variable
CEXeq1.3=1.3
CEXeq1.1=1.1
CEXeq1.0=1
CEXeq0.8=0.8
SPACEeq0.06=0.06
if (!is.null(output.fig.pdffile)) { pdf(output.fig.pdffile, width=6.5, height=sqrt(nrow(left.panel)), family=family.font) }
else {par(family=family.font)}
if (is.null(subtype.panel)) {
layout(mat=matrix(c(1, 2, 3, 4, 5, 6), nrow=2, byrow=TRUE), widths=c(40, 120, 40), heights=c(40, 80))
} else {
if (is.null(subtype.2.panel)) {
layout(mat=matrix(c(1, 2, 3, 4, 5, 6, 7, 8, 9), nrow=3, byrow=TRUE), widths=c(40, 120, 40), heights=c(40, 5, 80))
} else {
layout(mat=matrix(c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12), nrow=4, byrow=TRUE), widths=c(40, 120, 40), heights=c(40, 3, 3, 80))
}
}
### Set colors for legend, e.g. Purples, BuGn, BuPu
if (is.null(syn.nonsyn.legend.color)) {
syn.nonsyn.legend.color <- brewer.pal(3, "Purples")[2:3]
} else {
syn.nonsyn.legend.color <- syn.nonsyn.legend.color
}
if (sort.by.qvalue) {
right.panel <- right.panel[order(right.panel[,2], decreasing=F),]
genes <- right.panel[,1]
} else { # Sorting by mutation frequency
left.panel <- left.panel[order(left.panel[,2], decreasing=T),]
genes <- left.panel[,1]
}
genes <- as.character(genes)
middle.panel <- middle.panel[genes,]
rownames(left.panel) <- as.character(left.panel[,1])
rownames(right.panel) <- as.character(right.panel[,1])
right.panel <- right.panel[genes,]
left.panel <- left.panel[genes,]
middle.panel <- t(middle.panel)
# A dynamic way to order all of a matrix column, refer to following web link:
# http://stackoverflow.com/questions/12077413/order-a-matrix-by-multiple-column-in-r
middle.panel <- middle.panel[do.call(order, as.data.frame(-middle.panel)),]
middle.panel <- middle.panel[,rev(genes)]
sample.order <- rownames(middle.panel)
subtype.categs <- sort(unique(subtype.panel[,2]))
subtype.categs.num <- length(subtype.categs)
if (!is.null(subtype.panel)) {
rownames(subtype.panel) <- subtype.panel[,1]
subtype.data <- as.integer(subtype.panel[,2])
names(subtype.data) <- as.character(subtype.panel[,1])
subtype.data <- subtype.data[sample.order]
}
if (!is.null(subtype.2.panel)) {
rownames(subtype.2.panel) <- subtype.2.panel[,1]
subtype.2.data <- as.integer(subtype.2.panel[,2])
names(subtype.2.data) <- as.character(subtype.2.panel[,1])
subtype.2.data <- subtype.2.data[sample.order]
}
if (group.by.subtype) {
out <- list()
for (i in 1:subtype.categs.num) {
subtype.tmp.mat <- middle.panel[names(subtype.data[subtype.data==subtype.categs[i]]),]
subtype.tmp.mat <- subtype.tmp.mat[do.call(order, as.data.frame(-subtype.tmp.mat)),]
out[[i]] <- subtype.tmp.mat
}
middle.panel <- do.call('rbind', out)
sample.order <- rownames(middle.panel)
subtype.data <- subtype.data[sample.order]
if (!is.null(subtype.2.panel)) {
subtype.2.data <- subtype.2.data[sample.order]
}
}
### Plot an empty figure and legend for Syn. and Non-syn. mutations
plot.empty.graph <- function() plot(0,xaxt='n',yaxt='n',bty='n',pch='',ylab='',xlab='')
par(mar=c(0, 0, 0, 0))
plot.empty.graph()
#legend(x=legend.syn.nonsyn.x,y=legend.syn.nonsyn.y, legend=c("Syn.", "Non syn."), fill=syn.nonsyn.legend.color, box.col=NA, cex=1.3, border=syn.nonsyn.legend.color, x.intersp=0.2)
if (is.null(syn.nonsyn.legend.ncol)) {
legend('bottom', legend=c("Syn.", "Non syn."), fill=syn.nonsyn.legend.color, box.col=NA, cex=CEXeq1.3, border=syn.nonsyn.legend.color, x.intersp=0.2)
}
# Configure and plot for the upper panel
#par(mar=c(0, 0.16, 3, 0))
par(mar=upper.panel.mar)
rownames(upper.panel) <- as.character(upper.panel[,1])
upper.panel <- upper.panel[sample.order,]
rownames(upper.panel) <- NULL # Remove row names to prevent automatic barplot labeling
#par(mgp=c(0,-0.4,-1.3))
upper.panel.dat = t(upper.panel[,c(3,2)]/(32 * 1))
barplot(upper.panel.dat, border=NA, col=rev(syn.nonsyn.legend.color), las=2, axes=F, space=SPACEeq0.06, xaxs='i')
# Use tcl to control tick length, lwd to set axis width, col to set axis color
axis(side=2, cex.axis=CEXeq1.0, las=2, tcl=-0.3, lwd=0.6, col='grey20')
mtext(text="Mutations/Mb", side=2, padj=-3.5, cex=0.75)
if (!is.null(syn.nonsyn.legend.ncol)) {
if (syn.nonsyn.legend.ncol == 1) {
legend(syn.nonsyn.legend.loc, legend=c("Syn.", "Non syn."), fill=syn.nonsyn.legend.color, box.col=NA, cex=CEXeq1.3, border=syn.nonsyn.legend.color, x.intersp=0.2, horiz=F)
}
else {
legend(syn.nonsyn.legend.loc, legend=c("Syn.", "Non syn."), fill=syn.nonsyn.legend.color, box.col=NA, cex=CEXeq1.3, border=syn.nonsyn.legend.color, x.intersp=0.2, horiz=T)
}
}
if (is.null(image.color)) {
#image.color = brewer.pal(7,"Set1")
#image.color[6] = brewer.pal(7, "PuRd")[3]
#image.color[1] = brewer.pal(7, "BuPu")[1]
data('LandscapeColor', package="hustlxc")
image.color=image.color.lixc
}
#par(mar=c(0, 1, 2, 1))
par(mar=legend.panel.mar)
plot.empty.graph()
#legend(x=legend.panel.x, y=legend.panel.y, legend=legend.panel[,2][1:6], fill=rev(image.color), box.col=NA, cex=1.1, border=rev(image.color))
if (is.null(legend.panel.ncol)) {
legend("bottomleft", legend=legend.panel[,2][1:6], fill=rev(image.color), box.col=NA, cex=CEXeq1.0, border=rev(image.color))
}
if (!is.null(subtype.panel)) {
subtype.panel.color = subtype.panel.color[1:subtype.categs.num]
par(mar=c(0, 1, 0, 3))
plot.empty.graph()
par(mar=c(0, 1.3, 0., 1.16))
image(as.matrix(subtype.data), col=subtype.panel.color, axes=F)
grid(dim(middle.panel)[1], 1, col="white", lty=1, lwd=SPACEeq0.06)
mtext(side=2, text=subtype.panel.identifier, las=2, cex=subtype.panel.identifier.cex)
par(mar=c(0, 0, 0.0, 1))
plot(0, type='n', xlim=c(0,10), ylim=c(0, 1), axes=F)
#if (subtype.categs.num == 3) {
# text.width=c(0,3.5, 2.5)
#} else if (subtype.categs.num == 2) {
# text.width=c(0,3)
#}
# Use text.width to adjust distances between adjacent strings in the legend
# http://stackoverflow.com/questions/21262472/adjust-spacing-between-text-in-horizontal-legend
# Use inset to adjust distance between legend and margins, see ?legend
legend("left", text.width=subtype.legend.text.width, legend=subtype.labels, fill=subtype.panel.color, border=NA, box.col=NA, x.intersp=0.05, ncol=subtype.categs.num, xpd=T, inset=-0.04, cex=subtype.label.cex)
}
if (!is.null(subtype.2.panel)) {
subtype.2.categs <- sort(unique(subtype.2.panel[,2]))
subtype.2.categs.num <- length(subtype.2.categs)
subtype.2.panel.color = subtype.2.panel.color[1:subtype.2.categs.num]
par(mar=c(0, 1, 0, 3))
plot.empty.graph()
par(mar=c(0, 1.3, 0., 1.16))
image(as.matrix(subtype.2.data), col=subtype.2.panel.color, axes=F)
grid(dim(middle.panel)[1], 1, col="white", lty=1, lwd=SPACEeq0.06)
mtext(side=2, text=subtype.2.panel.identifier, las=2, cex=subtype.2.panel.identifier.cex)
par(mar=c(0, 0, 0.0, 1))
plot(0, type='n', xlim=c(0,10), ylim=c(0, 1), axes=F)
#if (subtype.2.categs.num == 3) {
# text.width=c(0,3,3)
#} else if (subtype.2.categs.num == 2) {
# text.width=c(0,3)
#}
# Use text.width to adjust distances between adjacent strings in the legend
# http://stackoverflow.com/questions/21262472/adjust-spacing-between-text-in-horizontal-legend
# Use inset to adjust distance between legend and margins, see ?legend
legend("left", text.width=subtype.2.legend.text.width, legend=subtype.2.labels, fill=subtype.2.panel.color, border=NA, box.col=NA, x.intersp=0.05, ncol=subtype.2.categs.num, xpd=T, inset=-0.04, cex=subtype.2.label.cex)
}
# The left panel
rownames(left.panel) <- as.character(left.panel[,1])
left.panel[,1] <- NULL
left.panel <- left.panel[rev(genes),]
#par(mar=c(5.5, 1.5, 0, 2.9))
par(mar=left.panel.mar)
#par(mgp=c(1, 0.5, 0))
if (is.null(left.panel.xlim)) {
left.panel.xlim = c(-1 * max(left.panel[,1]) - 1, 0)
}
bp <- barplot(t(-1 * left.panel[,1]), horiz=T, border=NA, col=rev(syn.nonsyn.legend.color), space=SPACEeq0.06, axes=F, xlim=left.panel.xlim, yaxs='i')
at <- seq(left.panel.xlim[1], 0, length.out=3)
#axis(side=1, at=c(0, -10, -20, -30), labels=c("0", "10", "20", "30"), cex.axis=1, xpd=TRUE)
axis(side=1, at=at, labels=as.character(at * -1), cex.axis=1, xpd=TRUE, lwd=0.6, tcl=-0.3, col="grey20")
#axis(side=1, cex.axis=1, xpd=TRUE)
axis(side=4, at=bp, labels=left.panel[,2], las=2, tick=F, hadj=0.01, cex.axis=left.panel.cex.axis)
mtext(text=left.panel.xlab, side=1, padj=1.6, font=1, cex=left.panel.xlab.cex)
#text(x=-60, y=8, labels="73", col="white")
# The middle panel
#par(mar=c(6, 1.3, 0.5, 1.16))
par(mar=middle.panel.mar)
image(middle.panel, col=image.color, axes=F)
grid(dim(middle.panel)[1], dim(middle.panel)[2], col="white", lty=1, lwd=SPACEeq0.06)
if (!is.null(legend.panel.ncol) && legend.panel.ncol==3) {
legend(x=0.1, y=legend.panel.y, legend=legend.panel[,2][1:6], fill=rev(image.color), box.col=NA, cex=CEXeq1.1, border=rev(image.color), xpd=T, ncol=3, x.intersp=0.2)
}
else if (!is.null(legend.panel.ncol) && legend.panel.ncol==1) {
legend(x=0, y=legend.panel.y, legend=legend.panel[,2][1:6], fill=rev(image.color), box.col=NA, cex=0.9, border=rev(image.color), xpd=T, horiz=T, x.intersp=0.2)
} else if (!is.null(legend.panel.ncol)){
legend(x=0.2, y=legend.panel.y, legend=legend.panel[,2][1:6], fill=rev(image.color), box.col=NA, cex=CEXeq1.1, border=rev(image.color), xpd=T, ncol=2, x.intersp=0.2)
}
# The right panel
q.value <- right.panel[,2]
smallest.q.value.clipped=FALSE
if (any(q.value == 0)) {
q.value[q.value==0] <- min(q.value) / 2
smallest.q.value.clipped=TRUE
} else {
q.value.n = length(q.value)
I=which(q.value == min(q.value))
second.min.q.value=sort(q.value)[2]/100
if (second.min.q.value/min(q.value) >100) {
q.value[I]=sort(q.value)[2]/100
smallest.q.value.clipped=TRUE
}
}
right.panel[,2] <- q.value
right.panel$qscore <- -log10(right.panel[,2])
if (smallest.q.value.clipped) {
warning('q-value was clipped to make visualization better!')
}
par(mgp=c(1, 0.5, 0))
par(mar=right.panel.mar)
right.panel.color <- right.panel.color
bp <- barplot(rev(right.panel$qscore), horiz=TRUE, col=right.panel.color, border=NA, space=SPACEeq0.06, axes=F, yaxs='i')
axis(side=2, at=bp, labels=rev(right.panel[,1]), las=2, font=3, tick=F, hadj=0.9, cex.axis=right.panel.cex.axis, lwd=CEXeq1.0)
#axis(side=1, at=c(0,4,8), labels=c("0",'4','8'), xpd=TRUE)
at <- ceiling(seq(0, max(right.panel$qscore), length.out=3))
axis(side=1, at=at, labels=as.character(at), xpd=TRUE, lwd=0.6, tcl=-0.3, col="grey20")
mtext(text=right.panel.xlab, side=1, padj=1.6, font=1, cex=right.panel.xlab.cex)
if (!is.null(right.panel.qscore.cutoff) && right.panel.qscore.lwd!=0) {
abline(v=right.panel.qscore.cutoff, lty=right.panel.qscore.lty, col=right.panel.qscore.col, lwd=right.panel.qscore.lwd)
}
if (!is.null(output.fig.pdffile)) {
dev.off()
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.