Tuesday, August 18, 2009

coverage

...
[maq-steps]
...
maq pileup -p [your bfa] [your map] > pileup.out

cut -f 2-4 pileup.out > croppedpileup.out

#then launch R
R
#following are R commands
data <-read.table(file="croppedpileup.out",sep="\t",header=F)
colnames(data)<-c("pos","consensus","coverage")
depth<-mean(data[,"coverage"])
# depth now has the mean (overall)coverage
#set the bin-size
window<-101
rangefrom<-0
rangeto<-length(data[,"pos"])
data.smoothed<-runmed(data[,"coverage"],k=window)
png(file="cov_out.png",width=1900,height=1000)
plot(x=data[rangefrom:rangeto,"pos"],y=data.smoothed[rangefrom:rangeto],pch=".", cex=1,xlab="bp position",ylab="depth",type="l")
dev.off()

No comments: