Difference between genomecov & coverage

[wd...@gate.sinica.edu.tw]

I checked "genomecov coverage" and seems no similar question yet.

This is what I tried:

1. use coverage to get base-pair-wise read coverage in sample.bam (attached) against a fake GFF3 file with only one object representing an entire chromosome

MyPrompt$ /home/wdlin/Tools/bedtools2/bin/bedtools coverage -d -s -split -abam sample.bam -b Ppatens_152_chromosome.gff3 | ./coverage.pl new

2. use genomecov to get base-pair-wise read coverage in sample.bam of the same chromosome

MyPrompt$ /home/wdlin/Tools/bedtools2/bin/bedtools genomecov -ibam sample.bam -d -split -strand + | perl -ne 'chomp; @token=split; print "$token[-1]\n"' > old.scaffold_1.coverage.txt

Seems some minor difference.

MyPrompt$ diff new.scaffold_1.coverage.txt old.scaffold_1.coverage.txt
1028201c1028201
< 6
---
> 5

I checked descriptions of coverage and genomecov, about -split they all said

        -split          Treat "split" BAM or BED12 entries as distinct BED intervals.
                        when computing coverage.
                        For BAM files, this uses the CIGAR "N" and "D" operations
                        to infer the blocks for computing coverage.
                        For BED12 files, this uses the BlockCount, BlockStarts, and BlockEnds
                        fields (i.e., columns 10,11,12).

and the read mapping is


Is this a bug? or actually genomecov and coverage compute depth differently?

Thank you in advanced.

Wendar

[Aaron Quinlan]

Hi Wendar,

As far as I know, they should indeed be identical if -split is used.  The one thing I notice is that your call to genomecov uses the “-strand +” option, which would restrict the count to solely those reads on the forward strand.  Given that your coverage from genomecov was less than that from coverage, could this explain the difference (I can’t tell by eye from the IGV shot you sent)?

Aaron

--
You received this message because you are subscribed to the Google Groups "bedtools-discuss" group.
To unsubscribe from this group and stop receiving emails from it, send an email to bedtools-discu...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
<sample.bam><coverage.pl><Ppatens_152_chromosome.gff3>

[wd...@gate.sinica.edu.tw]

No. The -s option of coverage is for "same strandedness" (-S for "different strandedness"), and the strand of the fake object is +. I checked all alignment in sample.bam, all reads inside are forward mapped.

Wendar