SEQanswers

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

Reply
 
Thread Tools
Old 09-27-2012, 06:15 AM   #1
krueml
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 - | 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
__________________
Thank you!

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

Join Date: Mar 2010
Posts: 56
Default

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
krueml
Junior Member
 
Location: US

Join Date: Aug 2012
Posts: 3
Default

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.

Code:
##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!

K
krueml 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 02:36 AM.


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