Plotting coverage from bamfiles

Under development.

Using samtools depth and R to plot the coverage of mapping from bamfiles. -aa option to include all nucleotides in reference, and -d to increase the maximum coverage depth.

# Using samtools v.1.3.1
# Setting max depth to 1000000 with -d (8000 default).
samtools depth -aa -d 1000000 input.bam | grep "contig_youwant_to_count" | gzip > coverage.txt.gz

# Enter R
library(tidyverse)
# Load data
cov <- read_tsv(file = "coverage.txt.gz", col_names = FALSE) %>% 
	rename("Position" = X2) %>% rename("Coverage" = X3)
# Plot
cov %>% select(Position, Coverage) %>% 
	ggplot(aes(Position, Coverage)) + 
	geom_line()

Leave a Comment