Script to calculate coverage of target regions for exome sequencing

coverage_stats.r

library(RColorBrewer)
rm(list=ls())
samples = read.table("Samples.txt",header=T,stringsAsFactors=F)
coverage = list()
cumulative_coverage = list()
for (i in 1:length(samples$files)) {
  coverage[[i]] <- read.table(samples$files[i])
  cumulative_coverage[[i]] <- 1-cumsum(coverage[[i]][,5])
}


cols <- brewer.pal(length(coverage), "Dark2")

pdf("coveragePlot.pdf",width=10,height=8)
plot(coverage[[1]][1:401,2], cumulative_coverage[[1]][1:401], type='n', yaxt="n", xlab="Depth", ylab="Fraction of capture target bases >= depth", ylim=c(0,1.0), main="Target Region Coverage (Morex Assembly)")
abline(v = 10, col = "gray60")
abline(v = 30, col = "gray60")
abline(v = 50, col = "gray60")
abline(v = 80, col = "gray60")
abline(v = 100, col = "gray60")
abline(h = 0.25, col = "gray60")
abline(h = 0.50, col = "gray60")
abline(h = 0.75, col="gray60")
abline(h = 0.90, col = "gray60")
axis(1, at=c(10,30,50,80), labels=c(10,30,50,80))
axis(2, at=c(0.90), labels=c(0.90))
axis(2, at=c(0.75), labels=c(0.75))
axis(2, at=c(0.50), labels=c(0.50))
axis(2, at=c(0.25), labels=c(0.25))

for (i in 1:length(coverage)) points(coverage[[i]][1:401, 2], cumulative_coverage[[i]][1:401], type='l', lwd=3, col = cols[i])
legend("topright", legend=samples$name, lty=1, lwd=4, col=cols)
dev.off()