This is the same problem I recently had, but I was not able to find a good solution. I was finally able to work a method which I am sharing here. While this question is 2 years old but I am hoping there are many individuals trying to work out the solution to this problem, so I hope this will be helpful for others.

Identification of the depth of coverage is quite useful in 1) identifying the regions that might have potential paralogous alignment 2) finding the coverage at your regions of interest.

Here is a step by step guide to the problem/question.

Step #1) First identify the depth at each locus from a bam file.

I have found samtools depth option more useful in this regard, when coverage at each locus is desired.

samtools depth deduped_MA605.bam > deduped_MA605.coverage

The output file 'deduped_MA605.coverage' file will have 3 columns (Chr#, position and depth at that position) like below.

Chr  position depth (this header will be absent though)
1       3980        66
1       3981        68
1      3982     67
1      3983        67
1      3984     68

Step #2) Now, select the coverage (depth) by locus for each chromosome and/or regions

We can use the coverage file to plot it in R. But, the file is so large that it will suck up almost all the memory. It better to split the coverage by chromosome (or region of the chromosome if required).

To select the coverage for a particular chromosome (Chr#1 in my case)

awk '$1 == 1 {print $0}' deduped_MA605.coverage > chr1_MA605.coverage

To select coverage from chr #2

awk '$1 == 2 {print $0}' deduped_MA605.coverage > chr2_MA605.coverage

If the chrosomosome has string characters it can be adjusted as

awk '$1 == "chr2" {print $0}' deduped_MA605.coverage > chr2_MA605.coverage

Step #3) To plot the data in R this coverage file will need to be imported and the headers need to be added.

MA605.chr2 <- read.table("/media/everestial007/SeagateBackup4.0TB/New_Alignment_Set/05-B-deDupedReads-forScanIndel/chr2_MA605.coverage",
         header=FALSE, sep="\t", na.strings="NA", dec=".", strip.white=TRUE)`

Note: The header of the column are automatically set as V1, V2 and V3.

To rename the headers

library(reshape) # loads the library to rename the column names

MA605.chr2<-rename(MA605.chr2,c(V1="Chr", V2="locus", V3="depth")) # renames the header

Now, plot the coverage by depth:

plot(MA605.chr2$locus, MA605.chr2$depth)

to get the wider/cleaner view of the plot use

library(lattice, pos=10) xyplot(depth ~ locus, type="p", pch=16, auto.key=list(border=TRUE), par.settings=simpleTheme(pch=16), scales=list(x=list(relation='same'), y=list(relation='same')), data=MA605.chr2, main="depth by locus - Chr2 (SampleMA605)")

The output plot looks like this:

Note:

• We can see that there is a region around the middle and next to it that has a very high coverage. This likely suggests paralogous alignment.

• Also, there is a gap in alignment. The region with no alignment is near (at) centromere where no consensus sequence has been found.

• Reads may be later filtered by coverage.

Thanks,

Step #1. Calculate the depth along the genomic locus by Samtools.

samtools depth reads.sort.bam > reads.sort.coverage

Step #2. Plot the coverage by ggplot function in R. (The demo data only includes one chromosome)

setwd("D://R")  # working path
coverage=read.table("reads.sort.coverage", sep="\t", header=F)
install.packages('reshape')
library(reshape)
coverage=rename(coverage,c(V1="Chr", V2="locus", V3="depth")) # renames the header
ggplot(coverage, aes(x=locus, y=depth)) +
geom_point(colour="red", size=1, shape=20, alpha=1/3) +
scale_y_continuous(trans = scales::log10_trans(), breaks = scales::trans_breaks("log10", function(x) 10^x))

The Picard suite is your ticket.

EDIT:

Specifically, CollectAlignmentSummaryMetrics gives you your alignment percentage (and many other metrics) for a BAM file input.

If you have a interval table (something like a BED file) which describes genomic regions of interest, you can calculate depth of coverage statistics using CalculateHSSummaryMetrics. This tool will also give you detailed alignment metrics as well, so you can run it in place of CollectAlignmentSummaryMetrics.

More generally, there's also BedTools' genomeCov and GATK's DepthOfCoverage

Take a look at this blog post from Stephen Turner., which has a step by step guide to visualize coverage for targeted NGS (exome) experiments

I wrote http://lindenb.github.io/jvarkit/WGSCoveragePlotter.html

I doubt you'll find an out of the box plotting solution for this but if all you need is coverage information, you could try using samtools mpileup to generate read coverage at each position in the BAM and make a plot using R.

I have used properly-paired qc-passed reads for chromosome wise depth calculation and plotting. However, still I can see 0 value for depth at some places in plot, even though I have used samtools depth and not samtools depth -a. Why is it so? Should I remove those 0 values from my data to proceed further or keep them. I am not clear on that.

Also samtools mpileup and samtools depth gives same output.

Please guide.

Thanks!

Using this commands, you can calculate breadth of coverage and depth of coverage. But only for one genome in the file bam. Can anyone suggest how to calculate for multiple genomes in one bam file? So that for each genome there is output of the width and depth of coverage. One by one it takes a long time to calculate. For example, outputting results for multiple genomes:

breadth of coverage  depth of coverage
genome1 70   200
genome2 80   300
...
genomen 15   100

breadth of coverage

./samtools depth -a ./My_data/aln_sort.bam | awk '{c++; if($3>0) total+=1}END{print (total/c)*100}'

average depth

./samtools depth -a /home/sergey/samtools-1.9/My_data/aln_sort.bam | awk '{c++;s+=$3}END{print s/c}'