equivalent to genomeCoverageBed

[lmanchon]

--Hi,

is there an equivalent command to this one with megadepth:

from bedtools i can output coverage per base with:

genomeCoverageBed -d -i chr20.bed -g chr20_size.genome

thank you --

[ChristopherWilks]

Hi @lmanchon,

While Megadepth doesn't have that specific capability in terms of output (and it doesn't operation on BED files as the source of reads), if you have a BAM file to start with, you could run the --coverage calculation and then extract out the per-base read counts you're interested in, as that information is there just consolidated over ranges of bases with the same read counts.

You'd get output similar to this example (note coordinates are in BED format starting at 0):

chr1    14008   14083   1
chr1    14083   14121   0
chr1    14121   14196   1

where the bases in range 14008-14083 are covered by 1 read.

This does require post-processing the Megadepth output yourself, something like this: megadepth <bam.file> --coverage | <script to extract out per-base read counts>

Also, you may want to pre-filter the first to the regions of interest using samtools or bamtools.

For more info, see: https://github.com/ChristopherWilks/megadepth#coverage-over-the-whole-genome