SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Newbie Question: Calculating physical coverage from genome coverage tristanstoeber Illumina/Solexa 4 06-24-2013 10:53 AM
Comparing Samples with 20Mbp coverage to 40MBP coverage sgoswami RNA Sequencing 0 03-09-2012 07:41 AM
differences in coverage BEDtools - samtools mpileup KNS Bioinformatics 3 12-01-2011 09:55 PM
coverage depth, bedtools and RNA-Seq experiment PFS Bioinformatics 2 03-18-2011 07:27 AM
low 454 coverage combined with high solexa coverage strob Bioinformatics 7 10-07-2010 10:14 AM

Reply
 
Thread Tools
Old 09-19-2013, 11:25 AM   #1
efoss
Member
 
Location: Seattle

Join Date: Jul 2011
Posts: 98
Default bedtools coverage is inaccurate

I am using bedtools to calculate coverage in an RNA-seq bam file. I get screwy readings when it encounters a gap where there are no reads. Rather than filling it in with zeros, it takes a value that is somewhere near the last read and uses that value to fill across the gap. For example, I just looked in IGV at a suspicious area in a bam file. IGV shows readcounts of 11, 11, 10, 10, 10, 9, 9, 9, 9, 1 and then a bunch of zeros. bedtools, on the other hand, lists the same numbers except that when it comes to the final 1, it says 7, and then the very long list of zeros is instead a very long list of 7s. Can anyone suggest a solution? The code I used is pasted below. (The last little thing where I pipe results to grep is to output only lines with non-zero read counts. Perhaps that's where I'm screwing up, but I don't think so.)

Thanks.

Eric



genomeCoverageBed -d -split -ibam input.bam -g /home/efoss/gene_databases/dot_genome_files/g1k_v37.genome | grep -v -P "\S+\t\S+\t0" > output.txt

Last edited by efoss; 09-19-2013 at 11:48 AM.
efoss is offline   Reply With Quote
Old 09-19-2013, 10:13 PM   #2
BAMseek
Senior Member
 
Location: St. Louis, MO, USA

Join Date: Apr 2011
Posts: 124
Default

Hey Eric,

My guess would be that if a read is split (such as spanning a splice junction), then that read may count towards the coverage for the entire interval, even the skipped bases. This link seems to talk about a similar issue. Looks like the solution in the post is to split those reads and try coverage again.

(Edit: although looks like you did use the -split argument so maybe it is a different issue).

Justin

Last edited by BAMseek; 09-19-2013 at 10:16 PM.
BAMseek is offline   Reply With Quote
Old 09-19-2013, 11:30 PM   #3
kokonech
Curious Character
 
Location: Berlin, Germany

Join Date: Sep 2010
Posts: 13
Default

Hi!

Make sure the BAM file is sorted by coordinate and the splitted alignment is defined by N cigar code (true for tophat by default).

Also you might check some alternative ways to calculate coverage, for example using coverageBed and intersection with annotation file:
http://okko73313.blogspot.de/2013/04...hromosome.html
kokonech is offline   Reply With Quote
Old 09-20-2013, 12:44 PM   #4
efoss
Member
 
Location: Seattle

Join Date: Jul 2011
Posts: 98
Default

Hi kokonech,

Thanks for the suggestion. I finally figured it out - there was a version of bedtools that had this as a known bug. I updated to the latest version and the problem went away.

I bookmarked your site. It looks very helpful.

Best wishes,

Eric
efoss is offline   Reply With Quote
Old 09-23-2013, 01:02 AM   #5
SylvainL
Senior Member
 
Location: Geneva

Join Date: Feb 2012
Posts: 177
Default

Hi efoss,

can you tell us which version of bedtools you were using before?

thanks
SylvainL is offline   Reply With Quote
Old 09-23-2013, 11:03 AM   #6
efoss
Member
 
Location: Seattle

Join Date: Jul 2011
Posts: 98
Default

Hi SylvainL,

I'm afraid I don't know which version it was. I think it was installed on our system in 2011. The one I'm using now is 2.17.0.

Eric
efoss is offline   Reply With Quote
Old 05-20-2014, 09:44 AM   #7
Nino
Member
 
Location: New York City

Join Date: Mar 2013
Posts: 27
Default

Hey Guys,

I am getting this weird error for my hg19.fa genome file.

this is my command line argument:
genomeCoverageBed -d -i Sample_23430_B14_T.sorted.dedup.rg.resorted.realigned.bed -g /pbtech_mounts/fdlab_store003/fdlab/genomes/human/hg19/indexes/bwa/genome/hg19.fa 1> Sample_23430_B14_T.sorted.dedup.rg.resorted.realigned.coverage

and I get this error:
Less than the req'd two fields were encountered in the genome file (/pbtech_mounts/fdlab_store003/fdlab/genomes/human/hg19/indexes/bwa/genome/hg19.fa) at line 1. Exiting.

anyone know why, I can't seem to figure it out.

Thanks!,
Nino
Nino is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 08:40 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO