plot_cover.R

##########################
# Author: B. Anderson
# Date: 21 Jan 2019
# Modified: Oct 2020, May 2021 (cumulative graph)
# Description: plot bp coverage given an input coverage file (arg 1) generated by samtools
##########################

args <- commandArgs()
cover_file <- args[6]


# read in the coverage file

data <- read.table(cover_file, header=FALSE, sep='\t', stringsAsFactors=FALSE)

contig_names <- unique(data[,1])	# unique row names in column 1


# start a pdf and run ploting for each contig
pdf('contig_cover.pdf', paper='special', width=11, height=6)


# if more than one contig, generate a cumulative plot
if (length(contig_names) > 1) {
	colnames(data) <- c('Contig', 'Position', 'Depth')
	data$Fill <- 1
	data$Cum <- cumsum(data$Fill)

	# calculate basic stats
	maxcov <- max(data$Depth)
	mincov <- min(data$Depth)
	meancov <- mean(data$Depth)
	sdcov <- sd(data$Depth)	

	# plot the coverage
	plot(data$Cum, data$Depth, type='l', pty='m',
		xaxt='n', xaxs='i', yaxs='i', xlim=c(1, max(data$Cum)), xlab='Cumulative Position', ylim=c(0, maxcov*1.1), ylab='Depth',
		main=paste0('All contigs, Coverage\nMax= ', maxcov, ', Min= ', mincov, ', Mean= ', 
		round(meancov, 2), ', StdDev= ', round(sdcov, 2)))

	# create an x axis appropriate to the scale
	power <- floor(log10(max(data$Cum))) - 1
	interval <- 10^(power)
	axis(1,	at=c(1, seq(interval, round(max(data$Cum), -power), interval), max(data$Cum)))

}

# Generate a plot per contig
index <- 1
for (name in contig_names) {

	# slice out a contig and rename the columns
	slice <- data[data[,1]==name,]		
	colnames(slice) <- c('Contig', 'Position', 'Depth')

	# calculate basic stats
	maxcov <- max(slice$Depth)
	mincov <- min(slice$Depth)
	meancov <- mean(slice$Depth)
	sdcov <- sd(slice$Depth)	

	# plot the coverage
	plot(slice$Position, slice$Depth, type='l', pty='m',
		xaxt='n', xaxs='i', yaxs='i', xlim=c(1, max(slice$Position)), xlab='Position', ylim=c(0, maxcov*1.1), ylab='Depth',
		main=paste0('Contig ', index, ' Coverage\nMax= ', maxcov, ', Min= ', mincov, ', Mean= ', 
		round(meancov, 2), ', StdDev= ', round(sdcov, 2)))

	# create an x axis appropriate to the scale
	power <- floor(log10(max(slice$Position))) - 1
	interval <- 10^(power)
	axis(1,	at=c(1, seq(interval, round(max(slice$Position), -power), interval), max(slice$Position)))

	index <- index + 1
}

dev.off()