Comparing measures of mtDNA mutation load in single cells.

mutationLoadMeasures.R

library(grDevices)
library(gtools)

dat = read.delim("data/RTdata.txt",sep="\t",stringsAsFactors=FALSE)
dat$PNUM = as.numeric(gsub("P","",dat$Patient))
dat$ID = sprintf("P%02d_%04d",dat$PNUM,dat$Cell.number)
dat$PAT = sprintf("P%02d",dat$PNUM)

colfunc = colorRamp(c("blue","yellow","red"),space="Lab")
colfun = function(x, alpha=1.0) {
 rgbcol = colfunc(x)
 return(rgb(rgbcol[,1]/255,rgbcol[,2]/255,rgbcol[,3]/255,alpha = alpha))
}

cats = colnames(dat)[4:6]
RC = c(15,16,17)
names(RC) = c("A","B","C")
RCalpha = 0.5
RCcols = c(rgb(1,0,0,RCalpha),rgb(0,0,1,RCalpha),rgb(0,0,0,RCalpha))
names(RCcols) = c("A","B","C")

perms = permutations(n=3,r=3,v=1:3)
pdf("MutationLoadEstimates.pdf",width=9,height=7.5)
for(i in 1:length(perms[,1])){
 x=perms[i,1]; y=perms[i,2]; z=perms[i,3]

 xvals = pmax(0,pmin(100,dat[,cats[x]]))
 yvals = pmax(0,pmin(100,dat[,cats[y]]))
 zvals = pmax(0,pmin(100,dat[,cats[z]]))

 xlab = paste(gsub("\\."," ",cats[x]),"(%)")
 ylab = paste(gsub("\\."," ",cats[y]),"(%)")
 zlab = paste(gsub("\\."," ",cats[z]),"(%)")

 zcols = colfun(zvals/100, alpha=0.7)
 zcols_full = colfun(zvals/100, alpha=1)
 zcol = rgb(0,0,0,0.3)

 layout(matrix(1:2,ncol=2), width = c(6,1),height = c(1,1))
 plot(xvals,yvals,pch=RC[dat$RC.cat], col=zcols,xlab=xlab,ylab=ylab,cex = 1.0,cex.lab=1.5,cex.axis=1.5)
 legend("bottomright",names(RC),col=zcol,pch=RC,bg="white")

 legend_image <- as.raster(matrix(colfun((100:0)/100,alpha=1), ncol=1))
 op = par(mar=c(0,0,3,0))
  plot(c(0,2),c(0,1),type = 'n', axes = F,xlab = '', ylab = '', main = zlab,cex.main=0.75)
  text(x=1.5, y = seq(0,1,l=5), labels = seq(0,100,l=5))
  rasterImage(legend_image, 0, 0, 1,1)
 par(op)

 zcols2 = colfun(zvals/100, alpha=1)
 op = par(mfrow=c(2,3))
  for (pat in sort(unique(dat$PAT))){
   plot(xvals[dat$PAT==pat],yvals[dat$PAT==pat],pch=RC[dat$RC.cat[dat$PAT==pat]], col=zcols2[dat$PAT==pat],xlab=xlab,ylab=ylab,cex = 1.0,cex.lab=1.5,cex.axis=1.5,main=pat,xlim=c(0,100),ylim=c(0,100))
  }
 legend("bottomright",names(RC),col=zcol,pch=RC,bg="white")
 par(op)

 layout(matrix(1:4,nrow=2,byrow=TRUE), width = c(6,6),height = c(6,6))
 for(cat in c("A","B","C")){
  op = par(mar=c(4,5,2.5,1))
  plot(xvals[dat$RC.cat.==cat],yvals[dat$RC.cat.==cat],pch=16, col=zcols[dat$RC.cat.==cat],xlab=xlab,ylab=ylab,cex = 1.0,cex.lab=1.5,cex.axis=1.5,main=cat,xlim=c(0,100),ylim=c(0,100))
  par(op)
 }

 legend_image <- as.raster(matrix(colfun((100:0)/100,alpha=1), ncol=1))
 op = par(mar=c(0,15,3,0))
  plot(c(0,2),c(0,1),type = 'n', axes = F,xlab = '', ylab = '', main = zlab,cex.main=0.75)
  text(x=1.5, y = seq(0,1,l=5), labels = seq(0,100,l=5))
  rasterImage(legend_image, 0, 0, 1,1)
 par(op)

 layout(matrix(1:2,ncol=2), width = c(6,1),height = c(1,1))
 plot(xvals,yvals,pch=16, col=RCcols[dat$RC.cat],xlab=xlab,ylab=ylab,cex = 0.5 + 1.5*zvals/100,cex.lab=1.5,cex.axis=1.5)
 legvals = seq(100,0,length.out=5)
 legend("bottomright",legend=legvals,pch=16,pt.cex=0.5 +1.5*legvals/100,col=zcol,cex=0.85,bg="white",title=zlab)
 legend("topright",names(RCcols),col=RCcols,pch=16,bg="white")
}
dev.off()