This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Vcf2Fq Not Calling Heterozygous Sites

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?

• 4,796 views • link
updated by Biostar Biostar • written by Justin Justin

1 answer

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.

• 245 views • link
updated by Ram Ram • written by Devon Ryan Devon Ryan

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.

• 0 views • link
written by Justin Justin

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

• 117 views • link
written by Devon Ryan Devon Ryan

Log in to answer this question.

More posts like this