genomeCoverageBed for paired end reads

Question

[Vanilla]

Hi all,

I'm now trying to used genomeCoverageBed in Bedtools v2.25.0 to calculated the genome coverage from bam files with paired-end reads. (http://bedtools.readthedocs.org/en/latest/content/tools/genomecov.html) However, it returns an error as:

ERROR: Unrecognized parameter: -pc

Does anyone know in which version of bedtools the "-pc" parameter is available? Why does it disappear in the latest version? If it is no longer available in bedtools, can anyone suggest a method to calculate the genome coverage for paired-end sequencing reads? Thanks very much for any comments and suggestions.

Best,
Vanilla

[James Ashmore]

Download the latest commit from github, it should be available.

[Vanilla]

Thanks James! Yeah it does work, even though the parameter is still not listed in the help message.

However, I don't think the output I got is correct. For example, the coverage for chr1:0-984444 I got without -pc parameter is (which is correct for raw reads visualized in IGV):

chr1    0    56228    0
chr1    56228    56247    2
chr1    56247    984425    0
chr1    984425    984444    2

Thus with -pc set, the coverage changed to be:

chr1    0    56228    0
chr1    56228    56247    1
chr1    56247    137447    0
chr1    137447    137455    1
chr1    137455    137456    17
chr1    137456    137459    20
chr1    137459    137460    32
chr1    137460    137464    35
chr1    137464    984425    36
chr1    984425    984444    37

The command line I used is:

genomeCoverageBed -ibam merge_noDuplicate_sort.bam -bga -pc > merge_NoDuplicate.bg

Is there anything wrong with my command? Thanks very mcuh!

Answer 1

[Carlo Yague]

I don't know about the -pc parameter, but to get the coverage of paired-end reads, samtools depth is quite straightforward.

[Vanilla]

Thanks Carlo! I haven't used samtools for that. May try it.

[Vanilla]

I just tried but still didn't resolve the "paired-end" issue...Could you clarify the usage of samtools depth for pared-end reads? Thanks!

[Carlo Yague]

First you need to split paired-end reads for the forward and reverse strand according to this tutorial, then use samtools on each bam file.