SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Read Depth in vcf (samtools / bcftools) Marie_Noir Bioinformatics 1 04-17-2012 06:48 AM
Maximum read Depth in Samtools Anjali Bioinformatics 2 01-16-2012 02:15 AM
How to determine maximum read depth in Samtools VarFilter Muhua wang Bioinformatics 0 10-06-2011 03:09 PM
How to decide maximum read depth in Samtools VarFilter Muhua wang Illumina/Solexa 0 10-06-2011 03:02 PM
What is read depth -D100 means in samtools? ketan_bnf Bioinformatics 1 07-28-2011 06:54 PM

Reply
 
Thread Tools
Old 05-16-2012, 12:32 PM   #1
vplagnol
Member
 
Location: London, UK

Join Date: Sep 2008
Posts: 13
Default samtools -d option to limit read depth

Dear all,

I need some help to understand the -d option of samtools, which is meant to do:

- At a position, read maximally INT reads per input BAM. [250]


I am working on PCR amplified data and the extreme read depth ~ 2,000 is affecting the QC steps. Hence I would like to make sure that I do not load more than 250 reads at each position. However, whatever I do using:

samtools mpileup -d 200 -f ...

I still get the same VCF output:
11 26463579 . C T 9.52 . DP=1990;VDB=0.0000;AF1=0.5;AC1=1;DP4=960,0,998,0;MQ=60;FQ=12.3;PV4=1,0,1,1
With the DP4 flags indicating a read depth around 2,000, so the full dataset without any cap on read depth. It looks like samtools is completely ignoring the -d option. Am I missing something here? How do I force samtools to cap the read depth?

Suggestions welcome,

Thanks,

Vincent
vplagnol is offline   Reply With Quote
Old 05-16-2012, 03:27 PM   #2
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

I've never used the -d option but I think you are misinterpreting what it does.

I bet what it does it looks at the bam, and after it sees that 200 reads have all started at the exact same base, it skips all the rest of the reads that start at that base, and moves onto the next base.

But that would still allow for lots more that 200x coverage at any one base.

So try dropping that -d even further. The other thing you can try is downsampling. Picard is one program that can do that.
swbarnes2 is offline   Reply With Quote
Old 05-16-2012, 03:48 PM   #3
vplagnol
Member
 
Location: London, UK

Join Date: Sep 2008
Posts: 13
Default

Thank you for your answer. But I can drop -d to very low levels, nothing changes in the resulting vcf. Moreover, because of the PCR amplification all my reads should start at the same place so the read depth should be effectively capped even with your interpretation. It looks to me like this samtools -d option is broken.

Thanks for the Picard suggestion, I think I'll try that, but it would be easier if it could be implemented at the calling stage rather than creating yet another BAM file.
vplagnol is offline   Reply With Quote
Old 05-21-2012, 07:19 PM   #4
me_myself_andI
Member
 
Location: Singapore

Join Date: Nov 2010
Posts: 30
Default

I might be wrong, but I think -d will only affect the pileup output, but not the snp prediction.
me_myself_andI is offline   Reply With Quote
Old 07-10-2013, 02:25 PM   #5
fongchun
Member
 
Location: Vancouver, BC

Join Date: May 2011
Posts: 55
Default

Quote:
Originally Posted by vplagnol View Post
Thank you for your answer. But I can drop -d to very low levels, nothing changes in the resulting vcf. Moreover, because of the PCR amplification all my reads should start at the same place so the read depth should be effectively capped even with your interpretation. It looks to me like this samtools -d option is broken.

Thanks for the Picard suggestion, I think I'll try that, but it would be easier if it could be implemented at the calling stage rather than creating yet another BAM file.
Did you ever find a solution to your problem with the -d parameter? I believe I am dealing with very similar data as you where I have PCR amplified amplicons and targeted re-sequencing generating depth of like 10,000. We you trying to limit the depth that is reported?

I am running into inconsistency where the depth reported by samtools is very small compared to what is reported in IGV.
fongchun is offline   Reply With Quote
Old 07-10-2013, 03:09 PM   #6
fongchun
Member
 
Location: Vancouver, BC

Join Date: May 2011
Posts: 55
Default

For my situation, I am running into inconsistencies with what is reported in IGV and samtools mpileup.

Where in IGV, I see 9162 reads supporting a position. When I run samtools mpileup like this:

Code:
samtools mpileup -A -B -q 1 -f ref.fa in.bam > out.mpileup
This position is reported as only containing 562 reads. I figured setting the -d parameter to a higher number would just solve this:

Code:
samtools mpileup -A -B -q 1 -d 10000 -f ref.fa in.bam > out.mpileup
But instead of returning 9162, it only reports 2562. When I increased it to 50,000:

Code:
samtools mpileup -A -B -q 1 -d 50000 -f ref.fa in.bam > out.mpileup
This ended up reporting all 9162 reads at this positions. I am not clear as to how the -d parameter actual works because it is behaving oddly. Also, by default it says the max depth is 250, but when you run samtools without using -d it always says the max depth per bam is set to 8000 at a position. So again, not sure why there is this discrepancy...

Does anyone know?

Thanks,
fongchun is offline   Reply With Quote
Old 07-10-2013, 06:37 PM   #7
BAMseek
Senior Member
 
Location: St. Louis, MO, USA

Join Date: Apr 2011
Posts: 124
Default

I would suggest setting the -d flag very high (like 100 million) in order to get accurate coverage results. If your coverage at a position is above the threshold, then the reported depth may not be correct. It may report a number much less than the threshold even if the coverage is above the threshold. This depth flag was probably put in there to limit the amount of memory being used - otherwise the maximum memory needed is proportional to the deepest depth. Setting the -d flag very high will make sure you don't hit that cap, which is especially useful for amplicons since you expect many reads to start at the same position.

Justin
BAMseek is offline   Reply With Quote
Old 07-11-2013, 07:41 AM   #8
fongchun
Member
 
Location: Vancouver, BC

Join Date: May 2011
Posts: 55
Default

Quote:
Originally Posted by BAMseek View Post
I would suggest setting the -d flag very high (like 100 million) in order to get accurate coverage results. If your coverage at a position is above the threshold, then the reported depth may not be correct. It may report a number much less than the threshold even if the coverage is above the threshold. This depth flag was probably put in there to limit the amount of memory being used - otherwise the maximum memory needed is proportional to the deepest depth. Setting the -d flag very high will make sure you don't hit that cap, which is especially useful for amplicons since you expect many reads to start at the same position.

Justin
Thanks for that suggestion. Will implement.
fongchun 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 06:02 PM.


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