Skip to content

Instantly share code, notes, and snippets.

@CnrLwlss CnrLwlss/mtDNA.R
Last active Dec 19, 2015

Embed
What would you like to do?
Generating a circular plot of sequencing coverage along the mitochondrial genome, annotating gene locations.
#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")
@VivekTodur

This comment has been minimized.

Copy link

commented Apr 29, 2015

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.