Go Back   SEQanswers > Bioinformatics > Bioinformatics

Similar Threads
Thread Thread Starter Forum Replies Last Post
haploid variant calling in SAMtools HTrewby Bioinformatics 12 01-13-2015 07:14 AM
minimum depth variant calling samtools/gatk m_elena_bioinfo Bioinformatics 1 12-06-2011 08:31 AM
variant calling kjaja Bioinformatics 1 11-04-2011 07:16 AM
Samtools variant calling questions Chiel Bioinformatics 2 06-07-2011 09:10 AM
variant calling using samtools -v- bcftools ksc Bioinformatics 2 04-13-2011 06:44 AM

Thread Tools
Old 09-27-2012, 06:15 AM   #1
Junior Member
Location: US

Join Date: Aug 2012
Posts: 3
Default variant calling with SAMtools

Hi, I am trying to repeat a papers steps to detect somatic variants. My earlier steps were using bwa, picard sortsam and addedreadgroup, and gatk indelrealigner. now I'm stucking at the variant calling step where the paper says I should use SAMtools varfilters with 3 different options. I found out how to put 2 of them but the 3rd one gives me headache since last week:

Minimum consensus quality set to 20

My pipeline looks like this for now:
samtools mpileup -M 255 -ugf ref.fasta realigned.bam | bcftools view - | varFilter -Q 43 | awk '$6>=50'> var.flt.vcf

I also found out on that there should be a column called consensus quality or phred scaled likelihood that genotype is wrong, respectively, but I cannot find that column in my data.
samtools mpileup -M 255 -f ref.fasta realigned.bam > raw.pileup
my data looks like this. And there arent so many columns (even with option -OABs with does not change):

1 569924 T 2 .. GI
1 569925 G 2 .. HI
Thank you!

krueml is offline   Reply With Quote
Old 09-27-2012, 11:51 PM   #2
Location: Zurich

Join Date: Mar 2010
Posts: 56

Hi krueml,

If you want to filter the bam file for a minimum mapping quality have a look at the -q option of samtools view:

samtools view -bq 20 your.bam > out.bam

If you want to see a "List of Phred-scaled genotype likelihoods" you should have a look at the output of your pipeline before using the awk command. The awk command also dismisses some header lines starting with '##' e.g.:
PHP Code:
##FORMAT=<ID=GL,Number=3,Type=Float,Description="Likelihoods for RR,RA,AA genotypes (R=ref,A=alt)">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="# high-quality bases">
##FORMAT=<ID=SP,Number=1,Type=Integer,Description="Phred-scaled strand bias P-value">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="List of Phred-scaled genotype likelihoods"> 
Btw: Your link points to the description of samtools pileup, but you use mpileup. Afaik there are some differences between both for example mpileup does not provide a consensus...
dschika is offline   Reply With Quote
Old 09-28-2012, 05:36 AM   #3
Junior Member
Location: US

Join Date: Aug 2012
Posts: 3

Yes sorry I posted the wrong link. This what I have posted was my mpileup output.
Well sorry that I wasn't clear. I was just asking how I can set the minimum consensus quality. After the whole pipeline, there is a FQ listed in the info of the vcf file.

##INFO=<ID=FQ,Number=1,Type=Float,Description="Phred probability of all samples being the same">

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  bar
1     9778579   .       C       T       62      .       DP=17;VDB=0.0000;AF1=1;AC1=2;DP4=0,0,13,0;MQ=60;FQ=-66  GT:PL:GQ        1/1:95,39,0:75
1     9781412   .       C       G       214     .       DP=66;VDB=0.0768;AF1=1;AC1=2;DP4=0,0,37,15;MQ=55;FQ=-184        GT:PL:GQ        1/1:247,157,0:99
1     9782556   .       C       T       222     .       DP=368;VDB=0.0599;AF1=1;AC1=2;DP4=2,0,263,42;MQ=59;FQ=-282;PV4=1,0.27,0.36,1    GT:PL:GQ        1/1:255,255,0:99
I changed the script by now of the varFilter but I wanted to know if there is another way (w/o changing the script), because in the paper there was no changing at all and again: they used the same version of samtools, so I was just curious.
Thank you!

krueml is offline   Reply With Quote

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 07:18 PM.

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