Question: Six phred-scaled genotype likelihoods for biallelic site
0
gravatar for hugo.zeberg
6 months ago by
hugo.zeberg0 wrote:

I have generated a VCF file from a bam file using samtools mpileup with the arguments -d 8000 -B -q30 -Q30.

I don't really understand the resulting PL field. Here is an example:

2   167050053   .   T   C,<*>   0   .   DP=3;I16=0,0,2,1,0,0,108,3888,0,0,111,4107,0,0,67,1539;QS=0,1,0;VDB=0.243476;SGB=-0.511536;MQSB=1;MQ0F=0    PL  101,9,0,101,9,101

In the above line there are 6 genotype likelihoods.

Based on the VCF specification. Biallelic sites has the ordering: AA,AB,BB; for triallelic sites the ordering is: AA,AB,BB,AC,BC,CC.

To me, it seems as if samtool treats the site as triallelic with 6 values, i.e. ,<*> is treat as an allele. In all positions were the reads match the reference I have <> and not the the dot ( "." ) I am used to.

I have used samtools 1.9.

Can anyone explain this behaviour? All help would be highly appreciated.

snp vcf genotype likelihood • 238 views
ADD COMMENTlink modified 6 months ago by finswimmer11k • written 6 months ago by hugo.zeberg0

So I found this:

5.5 Representing unspecified alleles and REF-only blocks (gVCF)

"A symbolic alternate allele <*> is used to represent this unspecified alternate allele"

So it indeed a triallelic site, with a unspecified alternate allele. How do I remove this unspecified alternate allele?

ADD REPLYlink written 6 months ago by hugo.zeberg0
0
gravatar for finswimmer
6 months ago by
finswimmer11k
Germany
finswimmer11k wrote:

Hello hugo.zeberg ,

the output you are showing is the output of mpileup and not the final vcf file. For this you need to use bcftools call.

When using samtools/bcftools for variant calling this is a two step procedure. In the first step mpileup is used for finding potential position where a variant can be. In the seconde step, only these position are investigate in more detail whether there this is a "real" variant.

So my understanding of the mpileup output is, that at this position there could be a REF (T), a specific ALT (C) or any other ALT (<*>). Add the callstep like this to get the final vcf:

$ bcftools mpileup -Ou -d 8000 -B -q30 -Q30 input.bam|bcftools call -mv -Ov -o output.vcf

Since v1.9 samtools mpileupis deprecated and should be replaced by bcftools mpileup.

fin swimmer

ADD COMMENTlink written 6 months ago by finswimmer11k
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: 1493 users visited in the last hour