SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
basic question about read groups efoss Bioinformatics 2 10-19-2011 05:32 PM
A basic base pair (dna) question ardmore General 12 08-24-2011 10:10 AM
depth of coverage basic question madsaan Bioinformatics 0 03-24-2011 07:40 AM
tophat - basic question - mapped reads chrisbala Bioinformatics 0 04-06-2010 10:36 AM
tophat - basic question - mapped reads chrisbala Illumina/Solexa 3 02-26-2010 09:38 AM

Reply
 
Thread Tools
Old 01-25-2012, 01:43 AM   #1
maria_mari
Member
 
Location: Germany

Join Date: Jan 2012
Posts: 17
Default a basic question about coverage

Hi everybody,
I have a basic question in NGS area
how can we calculate sequencing coverage (5X, 20X ...) at selected regions of interest? and what does it exactly mean?
It is calculated after sequencing and based on fastq file or after mapping to the genome?
maria_mari is offline   Reply With Quote
Old 01-25-2012, 01:50 AM   #2
mamons
Member
 
Location: Stockholm

Join Date: Nov 2011
Posts: 10
Default

Hello,

you need to map the reads first to know from what region they (hopefully) comes from.

One easy way to look for coverage in regions is to design a .bed file (http://genome.ucsc.edu/FAQ/FAQformat.html#format1) with your regions of interest and compare them to the mapped result with bedtools coverageBed (http://code.google.com/p/bedtools/).
mamons is offline   Reply With Quote
Old 01-25-2012, 06:46 AM   #3
mgogol
Senior Member
 
Location: Kansas City

Join Date: Mar 2008
Posts: 197
Default

Once you have coverage in terms of read count (from coverageBed), to get coverage like 5x, you'll have to do

( read count * read length ) / length of area in question

So if you have 5 reads that are 50 bases long in a region that's 100 bases long, your coverage will be

(5 * 50) / 100 = 2.5x.

You could calculate a number for the whole genome by adding all the chromosome lengths, or you could do individual chromosomes or genes or windows throughout the genome or whatever is interesting.
mgogol is offline   Reply With Quote
Old 01-25-2012, 02:29 PM   #4
aggp11
Member
 
Location: Wisconsin

Join Date: Jun 2011
Posts: 87
Default

You could also use "samtools depth" to find the coverage at each position in your target region which like mamons said, you could create a bed file for that.

Then you can just take all the coverage counts and get a summary using R to get the mean, median etc. depth of coverage in the region of interest
aggp11 is offline   Reply With Quote
Old 01-29-2012, 12:22 PM   #5
MBender
Junior Member
 
Location: Saarbruecken, Germany

Join Date: Jan 2011
Posts: 9
Default

Hi,

has someone of you ever tried to calculate the coverage (using coverageBed) of a whole genome sequencing experiment (about 40x average coverage) on a relatively large amount of genomic features (like all refSeq genes). I tried to perform this task on a multicore processor and 16 GB Ram memory. After three days of calculation and a constant memory consumption of about 14 GB i stopped the process. I used the following command:

./coverageBed -abam <bam_file> -b <refSeq_exons.bed6> -hist >> histogram.txt

Is that normal? Some ideas?

@aggp11: Which version of samtools are you using? The command depth seems not to be present in my version
MBender is offline   Reply With Quote
Old 01-29-2012, 02:40 PM   #6
maria_mari
Member
 
Location: Germany

Join Date: Jan 2012
Posts: 17
Default

Hi

http://sourceforge.net/mailarchive/f...samtools-devel


* Added the `depth' command to samtools to compute the per-base depth with a
simpler interface. File `bam2depth.c', which implements this command, is the
recommended example on how to use the mpileup APIs.
maria_mari is offline   Reply With Quote
Old 01-30-2012, 06:44 AM   #7
aggp11
Member
 
Location: Wisconsin

Join Date: Jun 2011
Posts: 87
Default

@Mbender: I am using samtools version 0.1.18 .

@maria_maria & Mbender: with this latest version of samtools, you don't even have to worry about bam2depth. The following command is an example of how samtools depth works:

samtools depth -q 30 -b exons.bed exome.bam > test_q_20.coverage

Output:
chr1 14468 39
chr1 14469 39
chr1 14470 37
chr1 14471 39
chr1 14472 35
chr1 14473 34

Where the third column is the # of q30 or more reads at the given position.

Thanks,
Praful
aggp11 is offline   Reply With Quote
Old 01-30-2012, 04:12 PM   #8
MBender
Junior Member
 
Location: Saarbruecken, Germany

Join Date: Jan 2011
Posts: 9
Default

Many thanks.

Using samtools depth seems to calculate the coverage in given genomic regions in a feasible amount of time. By the way, the low performance of coverageBed when working on a large amount of genomic intervals is a known issue.

http://groups.google.com/group/bedto...c7a5f40b28b5ef

Best,

Matthias
MBender 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:16 AM.


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