Entering edit mode
9.5 years ago
bioguy24
▴
230
I am trying to modify this very useful post by Michael James Clark which he wrote a perl script to calculate average coverage in a bam.
coverage.pl
($num,$den)=(0,0);
while ($cov=<STDIN>) {
$num=$num+$cov;
$den++;
}
$cov=$num/$den;
print "Mean Coverage = $cov\n";
What I am trying to do is use that code to calculate average coverage in a bam for targeted regions in a bed file.
/path/to/samtools view -b in.bam <genomic region> | /path/to/samtools mpileup - | awk '{print $4}' | perl ~/coverage.pl
Thank you :).
If you looking for alternative, both
bedtools
and bedmap frombedOps
can do it standalone.Thank you :)
Here is a nice blog post describing cumulative coverage across bait regions using bedtools. Very elegant.
Thank you, I will try bedtools. Does samtools have a similar feature?