variant calling with SAMtools

[krueml]

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 - | vcfutils.pl varFilter -Q 43 | awk '$6>=50'> var.flt.vcf


I also found out on http://samtools.sourceforge.net/cns0.shtml 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.
After
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

[dschika]

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

##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...

[krueml]

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.