R: Bioconductor: converage per base

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)

## Original 
# 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)

## More biology
# What if reads have gaps with cigar string like : 20M2D20M
# Treat "D" as "N"
# "I" and "S" not be considered default 
dat <- coverage(grglist(bm,drop.D.ranges=T))