bedtools genomecov produces negative coverage.

Question

[Carlo Yague]

Hi, I'm trying to compute fragment coverage from paired-end (50x2) RNA-seq data using bedtools genomecov with the -pc option. It looks like it is working well in most cases, however, I noticed that regions where few/no reads are mapped are sometimes reported with negative coverage. How is it possible ? Am I doing something wrong ?

Thank you for your time.

Carlo

Code :

#with -pc option, I get 6 negative values before position 18

bedtools genomecov -ibam aligned_reverse.bam -d -split -pc > rev.cov
head -n 20 rev.cov | tail -n 10

spike_in_T7 11  0   
spike_in_T7 12  -1  
spike_in_T7 13  -1  
spike_in_T7 14  -1  
spike_in_T7 15  -1  
spike_in_T7 16  -1  
spike_in_T7 17  -1  
spike_in_T7 18  0   
spike_in_T7 19  0   
spike_in_T7 20  0

/

#without -pc option, I see a read starting at position 18...

bedtools genomecov -ibam aligned_reverse.bam -d -split > rev.cov
head -n 20 rev.cov | tail -n 10

spike_in_T7 11  0   
spike_in_T7 12  0   
spike_in_T7 13  0   
spike_in_T7 14  0   
spike_in_T7 15  0   
spike_in_T7 16  0   
spike_in_T7 17  0   
spike_in_T7 18  1   
spike_in_T7 19  1   
spike_in_T7 20  1

PS : removing the -split option didn't change the issue.

PS2 : the reads from the bam file are properly paired.

PS 3 : Indeed, there is no read mapped before position 18

samtools view aligned_reverse.bam | awk '{print $3,$4,$10}' # print only template, pos, and sequence
spike_in_T7 18 GGGGAGGCGAGAACACACCACAACGAAAACGAGCAAAACCCGGTACGCA
spike_in_T7 40 ACGAAAACGAGCAAAACCCGGTACGCAACACAAAAGCGAACAACGCGAA
spike_in_T7 45 AACGAGCAAAACCCGGTACGCAACACAAAAGCGAACAACGCGAAAAAGG
spike_in_T7 61 TACGCAACACAAAAGCAAACAACGCGAAAAAGGACACCGAAGCGGAAGC
spike_in_T7 61 TACGCAACACAAAAGCAAACAACGCGAAAAAGGACACCGAAGCGGAAGC
spike_in_T7 105 GAAGCAAAGACAACCAACAGAAAACAACCGCAAACAAACGGGACCAGAC
spike_in_T7 107 AGCAAAGACAACCAACAGAAAACAACCGCAAACAAACGGGACCAGACAA
spike_in_T7 116 AACCAACAGAAAACAACCGCAAACAAACGGGACCAGACAACGCACCAGC

[ATpoint]

Maybe before wasting time on chasing ghosts: what is your version of bedtools? I remember there was an issue with the -pc parameter but I do not remember which version it was exactly. Just retry it with the most recent one, maybe that already solves it.

[Carlo Yague]

V2.26, I think it is the latest.

[Santosh Anand]

what the regions 12:16 look like if you try to see reads by samtools view? Is there any?

[Carlo Yague]

There is none, I edited my post accordingly