Okazaki.R

# Problem with bedpe file format as generated by bedtools
# https://code.google.com/p/bedtools/issues/detail?id=152
#SRR364781 -> wt_sample
#SRR364782 -> wt_replicate
#SRR364783 -> pol32_sample
#SRR364784 -> pol32_replicate

#source("https://bioconductor.org/biocLite.R")
#biocLite("ShortRead")
#biocLite("org.Sc.sgd.db")
library(ShortRead)
library(org.Sc.sgd.db)
chrlen <- org.Sc.sgdCHRLENGTHS
names(chrlen)[names(chrlen)=="chrM"]="chrMito"

refseq=c("chrI","chrII","chrIII","chrIV","chrV","chrVI","chrVII","chrVIII","chrIX","chrX","chrXI","chrXII","chrXIII","chrXIV","chrXV","chrXVI","chrMito")
names(refseq)=c("ref|NC_001133|","ref|NC_001134|","ref|NC_001135|","ref|NC_001136|","ref|NC_001137|","ref|NC_001138|","ref|NC_001139|","ref|NC_001140|","ref|NC_001141|","ref|NC_001142|","ref|NC_001143|","ref|NC_001144|","ref|NC_001145|","ref|NC_001146|","ref|NC_001147|","ref|NC_001148|","ref|NC_001224|")

outdir="plots"
dir.create(outdir)

colw=rgb(15/255,164/255,210/255)
colc=rgb(238/255,128/255,34/255)

for(sname in c("SRR364781","SRR364782","SRR364783","SRR364784")){

fname=file.path("alignments",paste(sname,"_fixed.bam",sep=""))
bm=readGAlignmentPairs(fname)

readlen=rep(0,length(bm))
centre=rep(0,length(bm))

posind=as.logical(strand(first(bm))=="+")
negind=as.logical(strand(first(bm))=="-")

lpos=abs(end(last(bm))-start(first(bm)))
lneg=abs(end(first(bm))-start(last(bm)))

readlen[posind]=lpos[posind]
readlen[negind]=lneg[negind]

cpos=round((end(last(bm))+start(first(bm)))/2.0)
cneg=round((end(first(bm))+start(last(bm)))/2.0)

centre[posind]=cpos[posind]
centre[negind]=cneg[negind]

png(file.path(outdir,paste(sname,"ReadLengths.png",sep="_")))
plot(density(readlen[as.logical(bm@first@strand=="+")]),col=colw,type="l",xlab="Length of read",xlim=c(0,1000))
points(density(readlen[as.logical(bm@first@strand=="-")]),col=colc,type="l")
points(density(readlen),type="l")
dev.off()

# Can check the size of all objects below (expect them to be zero)
wrongstrand=bm[strand(first(bm))==strand(last(bm))]
wrongchrom=bm[seqnames(first(bm))!=seqnames(last(bm))]
c(length(wrongstrand),length(wrongchrom))
range(readlen)

bmpos=bm[strand(first(bm))=="+"]
bmneg=bm[strand(first(bm))=="-"]

covpos=coverage(bmpos)
covneg=coverage(bmneg)

tofile=TRUE

for(chr in bm@first@seqinfo@seqnames){
	prettychr=refseq[chr]
	print(paste("Plotting for",prettychr,chr))
	lenchr=bm@first@seqinfo@seqlengths[bm@first@seqinfo@seqnames==chr]

	posbool=(as.logical(bm@first@seqnames==chr)&as.logical(bm@first@strand=="+"))
	negbool=(as.logical(bm@first@seqnames==chr)&as.logical(bm@first@strand=="-"))
	centrepos=centre[posbool]
	lenpos=readlen[posbool]
	centreneg=centre[negbool]
	lenneg=readlen[negbool]

	if(tofile) png(file.path(outdir,paste(sname,"_",prettychr,"_align%02d.png",sep="")),height=1000,width=4000,pointsize=35)

	rng=1:lenchr
	plot(centrepos,lenpos,cex=0.2,pch=16,ylim=c(-1000,1000),xlim=c(0,lenchr+1),col=colw,
	xlab="Chromosome coordinate of read centre",ylab="Read length",main=paste(prettychr,chr),xaxt="n",xaxs="i",bty="l")
	points(centreneg,-lenneg,cex=0.2,pch=16,col=colc)
	axis(side=1,at=axTicks(1),labels=formatC(axTicks(1),format="d",big.mark=","))

	plot(NULL,type="n",xlab="Chromosome coordinate",ylab="Coverage",main=paste(prettychr,chr),ylim=c(-2000,2000),xlim=c(0,lenchr+1),xaxt="n",xaxs="i",bty="l")
	polygon(c(head(rng,1),rng,tail(rng,1),head(rng,1)),c(0,as.vector(covpos[[chr]])[rng],0,0),col=colw,border=colw)
	polygon(c(head(rng,1),rng,tail(rng,1),head(rng,1)),c(0,-1*as.vector(covneg[[chr]])[rng],0,0),col=colc,border=colc)
	axis(side=1,at=axTicks(1),labels=formatC(axTicks(1),format="d",big.mark=","))
	
	if(tofile) dev.off()
}

}

# Generate coverage plots from .bed files in supplementary materials
library(rtracklayer)
bdir="beds"
beds=c("GSM835650_wt_sample.bed","GSM835651_wt_replicate.bed","GSM835652_pol32_sample.bed","GSM835653_pol32_replicate.bed")
refseq=c(refseq,"plasmid")
chrlen[["plasmid"]]=8000
for(bed in beds){
  print(paste("Plotting for",bed))
  reads=import.bed(con=file.path(bdir,bed))
  cspos=reads[(strand(reads)=="+")]
  csneg=reads[(strand(reads)=="-")]
  covpos=coverage(cspos)
  covneg=coverage(csneg)
  
  for(chrno in 1:18){
    png(file.path(outdir,paste(bed,"_",chrno,".png",sep="")),height=1000,width=4000,pointsize=35)
    lenchr=chrlen[refseq[chrno]]
    rng=1:lenchr
    plot(NULL,type="n",xlab="Chromosome coordinate",ylab="Coverage",main=paste(bed,refseq[chrno]),ylim=c(-400,400),xlim=c(0,lenchr+1),xaxt="n",xaxs="i",bty="l")
    polygon(c(head(rng,1),rng,tail(rng,1),head(rng,1)),c(0,as.vector(covpos[[chrno]])[rng],0,0),col=colw,border=colw)
    polygon(c(head(rng,1),rng,tail(rng,1),head(rng,1)),c(0,-1*as.vector(covneg[[chrno]])[rng],0,0),col=colc,border=colc)
    axis(side=1,at=axTicks(1),labels=formatC(axTicks(1),format="d",big.mark=","))
    dev.off()
  }
}