samtools -d option to limit read depth

[vplagnol]

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

[swbarnes2]

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.

[vplagnol]

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.

[me_myself_andI]

I might be wrong, but I think -d will only affect the pileup output, but not the snp prediction.

[fongchun]

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]

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:

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:

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:

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,

[BAMseek]

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

[fongchun]

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.