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