Question: bedtools genomecov produces negative coverage.
2
gravatar for Carlo Yague
4.1 years ago by
Carlo Yague5.2k
Canada
Carlo Yague5.2k wrote:

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
ADD COMMENTlink modified 4.0 years ago • written 4.1 years ago by Carlo Yague5.2k
1

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.

ADD REPLYlink written 4.0 years ago by ATpoint42k

V2.26, I think it is the latest.

ADD REPLYlink written 4.0 years ago by Carlo Yague5.2k
1

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

ADD REPLYlink written 4.0 years ago by Santosh Anand5.2k

There is none, I edited my post accordingly

ADD REPLYlink written 4.0 years ago by Carlo Yague5.2k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1717 users visited in the last hour