Coverage.R

library(ShortRead)

# Root of .bam filename
sname="HG00100.unmapped.ILLUMINA.bwa.GBR.low_coverage.20130415"

# Construct path to file
fname=file.path("alignments",paste(sname,".bam",sep=""))

# Read alignments
bm=readGappedAlignments(fname)

# Calculate coverage at each position on each chromosome
dat=coverage(bm)

# Take the reads for chromosome 11 and thin (otherwise plots will take a lot of memory)
chrname="11"
rawreads=as.numeric(dat[[chrname]])
delt=10000
inds=seq(1,length(rawreads),delt)

# Prepare blank plot to show coverage along length of chromosome
chrid=paste("Chromosome",chrname)
plot(1:10,1:10,type="n",xaxt="n",xlim=c(0,134),ylim=c(0,120),ylab="Coverage",xlab="Chromosome Coordinate (Mb)",main=chrid)
abline(v=seq(0,1550,12.5),col="grey")
axis(1,at=seq(0,250,25))

# Scale coordinate number (from bases to Mb) and draw on plot
points(inds/1000000,rawreads[inds],type="l",lwd=1)

Good tutorial ~ p.s. use <- instead of = since this is R

Hi
I have trouble with these commands. Actually it works fine with your dataset which is chr11 .bam file.
But, when I use my dataset which is accepted_hits.bam from tophat output (it has all chr starting from 1 to y), it gives error (please check the steps 9 and 13):

1> library("ShortRead")
2> sname="accepted_hits"
3> fname=file.path("alignments",paste(sname,".bam",sep=""))
4> bm=readGAlignments(fname)
5> dat=coverage(bm)
6> chrname="11"
7> rawreads=as.numeric(dat[[chrname]])
8> delt=10000
9> inds=seq(1,length(rawreads),delt)
Error in seq.default(1, length(rawreads), delt) :
wrong sign in 'by' argument
{Then I changed 1 to 0.
9.1> inds=seq(0,length(rawreads),delt)
10> chrid=paste("Chromosome",chrname)
11> plot(1:10,1:10,type="n",xaxt="n",xlim=c(0,134),ylim=c(0,120),ylab="Coverage",xlab="Chromosome Coordinate (Mb)",main=chrid)
11> abline(v=seq(0,1550,12.5),col="grey")
12> axis(1,at=seq(0,250,25))
13> points(inds/1000000,rawreads[inds],type="l",lwd=1)
Error in xy.coords(x, y) : 'x' and 'y' lengths differ
????????
Can somebody help me please?