Question: Vcf2Fq Not Calling Heterozygous Sites
0
gravatar for Justin
5.5 years ago by
Justin440
United States
Justin440 wrote:

When I run the following command:

samtools mpileup -g -q 50 -Q 30 -r chr1:19950062-19950062 file.bam | bcftools view -c -

I get C/T as output in the 5th field:

#CHROM    POS    ID    REF    ALT    QUAL    FILTER    INFO    FORMAT
chr1    19950062    .    N    C,T    69    .    DP=59;QS=0.000000,0.678448,0.321552,0.000000;VDB=1.639595e-01;AF1=1;AC1=2;DP4=0,0,59,0;MQ=50;FQ=-129    PL    255,255,153,243,0,186

But when I run the same command but pipe it in vcfutils:

samtools mpileup -g -q 50 -Q 30 -r chr1:19950062-19950062 file.bam | bcftools view -c - | vcfutils.pl vcf2fq

The nucleotide at chr1:19950062 says C, and not C,T (or Y)

Is this a bug or am I doing something wrong?

• 2.8k views
ADD COMMENTlink modified 2.8 years ago by Biostar ♦♦ 20 • written 5.5 years ago by Justin440
1
gravatar for Devon Ryan
5.5 years ago by
Devon Ryan89k
Freiburg, Germany
Devon Ryan89k wrote:

It says C there because the FQ is negative. From the mpileup documentation:

If negative, FQ equals the minus phred-scaled probability of all chromosomes being identical. Notably, given one sample, FQ is positive at hets and negative at homs.

Edit: I should add that even if the FQ weren't negative, it still wouldn't do what you likely want. In those cases, it uses the REF and the first ALT allele to form the heterozygous call. This probably isn't ideal in the example that you posted.

ADD COMMENTlink modified 5.5 years ago • written 5.5 years ago by Devon Ryan89k

Oh I see, so if my REF was C and my ALT was T, it would output Y? I guess that means I should specify a reference fasta in samtools mpileup.

ADD REPLYlink written 5.5 years ago by Justin440
1

Correct, you'll get some not-so-useful results without supplying the fasta file.

ADD REPLYlink written 5.5 years ago by Devon Ryan89k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1390 users visited in the last hour