SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
samtools mpileup using -l or -r vplagnol Genomic Resequencing 3 05-08-2013 10:23 PM
samtools mpileup: -B or not to -B? yenhuahuang1 Bioinformatics 4 03-19-2013 03:21 AM
samtools mpileup Ling1 Bioinformatics 1 05-24-2012 10:09 AM
Many Ns in samtools mpileup mikheyev Bioinformatics 5 03-22-2012 09:56 AM
samtools mpileup hv4 Bioinformatics 2 10-29-2011 01:21 PM

Reply
 
Thread Tools
Old 06-13-2012, 03:08 AM   #1
bair
Member
 
Location: London

Join Date: Jan 2010
Posts: 65
Default samtools mpileup

Hello All,

I have a bam file generated by MiSeq + GATK (Realignment/recalibration), I use samtools (v0.1.18) / GATK for variant calling. For a site, I find a big difference about the depth between these two methods.

bcftools view sample.bcf | awk '$1=="chr7" && $2==150648198'
chr7 150648198 . A G,X 0 . DP=8;I16=1,0,0,1,24,576,25,625,60,3600,60,3600,25,625,25,625;VDB=0.0000 PLP:SP 19,0,18,22,21,40:2:0

Here DP=8 before quality control and only two reads passed quality control, thus samtools does not call variant for this site.

GATK with default setting calls this site as SNP with DP=516, and IGV shows there are 233A + 283G at this site with good quality (mapQ~37, baseQ~30.

This site is a true SNP confirmed by Sanger sequencing, so any explanation will be appreciated.

Many thanks!
bair is offline   Reply With Quote
Old 06-13-2012, 08:30 AM   #2
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

One quick thing to try; remake bcf, using -B in the samtools mpileup stage.

Sometimes, the BAQ calculations will count good variant reads as misalignments, and it will give them poor quality, which might explain your low coverage.
swbarnes2 is offline   Reply With Quote
Old 06-13-2012, 01:36 PM   #3
bair
Member
 
Location: London

Join Date: Jan 2010
Posts: 65
Default

Quote:
Originally Posted by swbarnes2 View Post
One quick thing to try; remake bcf, using -B in the samtools mpileup stage.

Sometimes, the BAQ calculations will count good variant reads as misalignments, and it will give them poor quality, which might explain your low coverage.
Thanks.

I remade bcf with -B, no change. DP=8.

Worked on the raw bam file without GATK realign/recalibration, samtools gave the same DP with/without -B option.
bair is offline   Reply With Quote
Old 08-30-2013, 10:22 AM   #4
jcgrenier
Member
 
Location: Montreal

Join Date: May 2011
Posts: 12
Default

Hi,

I had the same exact problem with my sample. I thought at the beginning that it was due to the soft clipping in my reads, but when I correct for that, I still have the same result. I have also played with the mapping and base qualities, still not able to recover the same depth for both programs.

Then, I tried the option -A in mpileup and pouff! Magic, I got the same numbers. I suppose that it's little bit risky to take those reads, but at least now I know what was doing mpileup.

Cheers
jcgrenier is offline   Reply With Quote
Old 11-21-2013, 09:33 AM   #5
petrie
Junior Member
 
Location: California

Join Date: Dec 2009
Posts: 5
Default

jcgrenier - could you tell me how you corrected for soft clipping? I have a lot of mysterious soft clipping in my bwa sw alignment, and I'd love to get rid of it. Just asking on the off chance that you had the same problem as me…Thanks!

See also: http://seqanswers.com/forums/showthr...119#post121119
petrie is offline   Reply With Quote
Old 11-21-2013, 11:07 AM   #6
jcgrenier
Member
 
Location: Montreal

Join Date: May 2011
Posts: 12
Default

Hi petrie,

For the moment, I'm kind of exploiting a bug in GATK to help me to solve this problem. If you go on their forum there :

http://gatkforums.broadinstitute.org...-option#latest

you will see why it does that. Anyways, here's the option that I used to kind of revert the softclipped bases. The problem was coming from the fact that the starting position of a softclipped read is in the sam file line not the one from the beginning of the read, but the first non-clipped base.

java -Djava.io.tmpdir=$SCRATCH -Xmx10g -jar ~/GenomeAnalysisTK-2.4-9-g532efad/GenomeAnalysisTK.jar -T ClipReads -l INFO -I $SCRATCH/$f.bam -o $f.hardclip.bam -R ~/References/hg19/ref.fasta -CR HARDCLIP_BASES -os $f.hardclip.stats

Cheers.
jcgrenier is offline   Reply With Quote
Old 11-21-2013, 11:11 AM   #7
jcgrenier
Member
 
Location: Montreal

Join Date: May 2011
Posts: 12
Default

Hi petrie,

For the moment, I'm kind of exploiting a bug in GATK to help me to solve this problem. If you go on their forum there :

http://gatkforums.broadinstitute.org...-option#latest

you will see why it does that. Anyways, here's the option that I used to kind of revert the softclipped bases. The problem was coming from the fact that the starting position of a softclipped read is in the sam file line not the one from the beginning of the read, but the first non-clipped base.

java -Djava.io.tmpdir=$SCRATCH -Xmx10g -jar ~/GenomeAnalysisTK-2.4-9-g532efad/GenomeAnalysisTK.jar -T ClipReads -l INFO -I $SCRATCH/$f.bam -o $f.hardclip.bam -R ~/References/hg19/ref.fasta -CR HARDCLIP_BASES -os $f.hardclip.stats

Cheers.
jcgrenier is offline   Reply With Quote
Old 11-23-2013, 03:10 PM   #8
petrie
Junior Member
 
Location: California

Join Date: Dec 2009
Posts: 5
Default

jcgrenier - that looks like it will solve my problem (either the hardclip_bases workaround or the revert_softclipped_bases if it's been updated). Thanks so much!
petrie is offline   Reply With Quote
Reply

Tags
false negative

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:59 AM.


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