Hello all,
I'm using samtools in order to get the depth of a file called file.bam, but the output I have is like this:
samtools depth file.bam > output.txt
11 88911833 2
11 88911834 2
11 88911835 2
I understand that the first column is the name of the reference sequence, the second column is the base index within the reference, and the third column is the depth of coverage for that base.
But I need at the end the read depth represented like 100X, 150X and something like this.
Is there any calculation should I follow/apply to the output file?
2 answers
You can calculate the average coverage (for covered bases):
samtools depth *bamfile* | awk '{sum+=$3} END { print "Average = ",sum/NR}'
This would be average coverage for all covered regions. If you want to include regions that were not covered in this calculation you need ot use: samtools depth -a
To get your average X coverage you would need to divide by the total size of your genome, instead of dividing by NR in the command above.
The total size can be calculated like this:
samtools view -H *bamfile* | grep -P '^@SQ' | cut -f 3 -d ':' | awk '{sum+=$1} END {print sum}'
Then you would redo the calculation substituting this in for NR.
I stole all of this from: Tools To Calculate Average Coverage For A Bam File?
You need to divide by the total number of bases you expect to be covering with your sequencing reads to get X.
It will work with any type of sequencing. By default if you follow the answer line by line, the sum will take into account all bases covered in your BAM file. If you did WES, this is likely going to include off-target sites so it will probably deflate your average coverage value. You need to manually replace NR (as described below in the code block) with the total number of bases covered by your WES kit or more ideally provide -b to the samtools depth command and give it a BED file of your WES target regions.
I downloaded a Regions.bed file from the Agilent Suredesign site and used it to estimate coverage by BEDtools. I didn't know about the "samtools depth -b BEDfile.bed" option so thank you very much! The coverage values in SAMtools with a BED file are roughly 9-10 points higher than the results from BEDtools (coverage around 50-60). So this is the actual average X coverage I believe. Again thanks a bunch!
Log in to answer this question.
-
bedcov output meaning
written by hello_bioinf •
-
Read depth per base DNA seq
written by Sabeen •
-
samtools depth -r
written by lz1903 •
-
Samtools depth to get the read depth
written by bandanaschapagain •
-
average depth across samples
written by ijlal.hyder2012 •
-
large bam file and get small mpileup file
written by bingnas •
-
Compare one of the columns of two seperate files and write the raw of the second file next to first…
written by emre •
-
[samtools depth] Calculate percent of the genome covered
written by madkitty
-
Produce A Bed Or Position List File Containing A List Of Regions For Mpileup
written by rob234king
-
Parsing Base String From Samtools Mpileup
written by Noushin N
I think you want the average depth of your sequencing data, to do that you can simply take the distribution of individual bases depths and compute the mean.