Generating a circular plot of sequencing coverage along the mitochondrial genome, annotating gene locations.

mtDNA.R

#source("http://bioconductor.org/biocLite.R")
#biocLite("ecolitk")
#biocLite("ShortRead")
library(ecolitk)
library(plotrix)
library(ShortRead)

# Read in some alignments
getReads=function(fname,rad=1.5,readmax=29000,logb=10){
	bm=readGappedAlignments(fname)
	dat=coverage(bm)
	rawreads=as.numeric(dat$mtDNA)
	kbmax=length(rawreads)
	metadata=strsplit(fname,"_")[[1]]

	lreads=log(pmax(1,rawreads),logb)
	ml=min(lreads)
	lreads=lreads-ml
	lreads=lreads/(log(readmax,logb)-ml)
	reads=rawreads/readmax

	pos=0:(kbmax-1)
	thets=pi/2+2*pi*pos/kbmax
	dxl=(lreads+rad)*cos(thets)
	dyl=(lreads+rad)*sin(thets)
	dx=(reads+rad)*cos(thets)
	dy=(reads+rad)*sin(thets)
	df=data.frame(Position=pos,Reads=rawreads,Scaled=reads,LogScaled=lreads,Theta=thets,X=dx,Y=dy,Xl=dxl,Yl=dyl,Barc=metadata[2])
	return(df)
}

# Draw a sunburst plot of mitochondrial genome
drawDNA=function(tracelst,clist,geneData,cexnum=0.5,cexlab=0.5,rado=1.5,radi=1.0,logged=FALSE,logb=10,readmax=29000,legends="",mlab="mtDNA",deltatheta=0,filled=TRUE){

	# Shift text by small angle to make it centred on the region it's labelling...
	offset=2*pi*deltatheta/360
	op=par(mai=c(0,0,0,0),bg = "white")
	# Draw a circle
	plot(NULL,xlim=c(-2.7,2.7),ylim=c(-2.7,2.7),axes=FALSE,ann=FALSE)


	# Trace scale
	rads=(10^seq(0,4))
	lrads=(log(rads,logb)-log(1,logb))/(log(readmax,logb)-log(1,logb))
	radslin=seq(0,readmax,500)
	if(!logged){rads=radslin; lrads=radslin/readmax}
	for (j in 1:length(rads)){
		rd=lrads[j]
		linesCircle(rado+rd,lty=2)
		text(0,1.025*rado+rd,rads[j],cex=0.5)
	}
	# Draw histogram
	tlist=c()
	for(i in 1:length(tracelst)){
		tr=tracelst[[i]]
		col=clist[i]
		tx=as.character(tr$Barc[1])
		tlist=c(tlist,tx)
		linew=2
		if(filled){
			if(logged){polygon(tr$Xl,tr$Yl,col=col,border=NA)}else{polygon(tr$X,tr$Y,border=NA,col=col)}
		}else{
			if(logged){points(tr$Xl,tr$Yl,type="l",col=col,lwd=linew)}else{points(tr$X,tr$Y,type="l",col=col,lwd=linew)}
		}
	}
	draw.circle(0,0,rado,nv=100,border="white",col="white",lty=1,lwd=1)
	
	for(i in 1:dim(geneData)[1]){
		if(geneData$Ending[i]<geneData$Starting[i]){
			end=geneData$Ending[i]+kbmax
		}else{
			end=geneData$Ending[i]
		}
		polygonArc(pi/2+2*pi*geneData$Starting[i]/kbmax,pi/2+2*pi*end/kbmax,radi,rado,edges=1000,col=geneData$col[i],border="white",lwd=1)
	}
	# Draw gene labels
	for(i in 1:dim(geneData)[1]){
		if(geneData$Ending[i]<geneData$Starting[i]){
			thetamid=pi/2+2*pi*((geneData$Starting[i]+kbmax+geneData$Ending[i])/2)/kbmax
		}else{
			thetamid=pi/2+2*pi*((geneData$Starting[i]+geneData$Ending[i])/2)/kbmax
		}
		thetaraw=pi/2+2*pi*geneData$Label[i]/kbmax
		if((thetaraw>=pi/2)&(thetaraw<3*pi/2)) {
			theta=thetaraw+offset
			ttext=pi+theta
			text(0.95*radi*cos(theta),0.95*radi*sin(theta),geneData$Shorthand[i],col="black",cex=cexlab,srt=360*ttext/(2*pi),adj=c(0,0))
		}else{
			theta=thetaraw-offset
			ttext=theta
			dist=0.95*radi-strwidth(as.character(geneData$Shorthand[i]))*cexlab
			text(dist*cos(theta),dist*sin(theta),geneData$Shorthand[i],col="black",cex=cexlab,srt=360*ttext/(2*pi),adj=c(0,0))
		}
		if(geneData$Line[i]=="Yes") {
		points(c(0.96*radi*cos(thetaraw),radi*cos(thetamid)),c(0.96*radi*sin(thetaraw),radi*sin(thetamid)),type="l",lwd=1)
		}
		
	}
	# Draw coordinate numbers
	kbs=seq(0,kbmax,1000)
	trad=1.025*radi
	for(kb in kbs){
		theta=pi/2+2*pi*(kb/kbmax)
		if((theta>=pi/2)&(theta<3*pi/2)) {
			theta=theta+offset
			ttext=pi+theta
			strw=strwidth(paste(kb/1000," kb",sep=""),cex=cexnum)
			dist=trad+strw
			text(dist*cos(theta),dist*sin(theta),paste(kb/1000," kb",sep=""),col="white",cex=1.0*cexnum,srt=360*ttext/(2*pi),adj=c(0,0))
		}else{
			ttext=theta
			text(trad*cos(theta),trad*sin(theta),paste(kb/1000," kb",sep=""),col="white",cex=1.0*cexnum,srt=360*ttext/(2*pi),adj=c(0,0))	
		}
	}
	text(0,0,mlab,cex=1.5)

	if(legends!="") {
		tlist=legends
		legend("topright",legend=tlist,lwd=3,col=clist,bg="white")
	}
	par(op)
}

# Read in alignment file and draw plot
mtDNAPlot=function(annotation,alignment,readmax=10000,histColour="blue",logged=TRUE,label=""){
	outerRad=1.5
	# Read in some mtDNA annotations
	mtDat=read.delim(annotation,sep="\t",header=TRUE,stringsAsFactors=FALSE)
	mtDat$col=mtDat$Colour
	N=dim(mtDat)[1]
	# Read in alignment files
	bamData=getReads(alignment,logb=10,rad=outerRad,readmax=readmax)
	# Make plot
	drawDNA(list(bamData),histColour,logged=logged,cexlab=.7,cexnum=0.675,rado=outerRad,radi=1.2,mlab=label,deltatheta=1.5,logb=10,readmax=readmax,geneData=mtDat)
}

mtDNAPlot("mtDNA_Colours.txt","IonXpress_018_R_2012_12_11_09_08_56_user_SN1-54_Auto_user_SN1-54_53.bam",readmax=1100,histColour="blue",logged=TRUE,label="CHiP-seq")

Hi, Thanks for the script. I am trying to execute, but I am getting error message saying, object 'kbmax' not found. Can you please help me to rectify...?

Thanks,
Vivek