Vcf2Fq Not Calling Heterozygous Sites

Question

[Justin]

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?

Answer 1

[Devon Ryan]

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.

[Justin]

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.

[Devon Ryan]

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