Question on calling SNPs using samtools/bcftools

[nkwuji]

Hi,

I am using samtools/bcftools to call SNPs. I have a question on bcftools, I don't understand why this SNP is not called by bcftools. For example, this is the result of test.bcf file, which is created by:

samtools mpileup -ugf ref.fa seq.bam >test.bcf

chr21 8588396 . T A,X 0 . DP=13;I16=1,9,0,3,378,14392,73,1809,600,36000,180,10800,93,1075,42,642 PL 27,0,57,190,199,232

From PL value, P(D|TT)=10^(-2.7), P(D|TA)=0. So TA allele should be the most possible allele.

But in the processed vcf using this command:

bcftools view -c test.bcf > test_count.vcf

This SNP is:

chr21 8588396 . T . 4.76 . DP=13;AF1=0.4995;CI95=0.5,0.5;DP4=1,9,0,3;MQ=60;PV4=1,5.9e-05,1,1 PL 27

There is no alternative base information, but the AF1 is about 0.5, which means this is an heterozygous allele. And, if I use the following command, there is no SNP called at all.

bcftools view -vcg test.bcf > test_count.vcf

Why the alternative base is not A??? Is there any option I can choose to call this SNP?

BTW, I am using samtools version 0.1.10-6.

Cheers,

Jun

[nkwuji]

Originally Posted by nkwuji

Hi,

I am using samtools/bcftools to call SNPs. I have a question on bcftools, I don't understand why this SNP is not called by bcftools. For example, this is the result of test.bcf file, which is created by:

samtools mpileup -ugf ref.fa seq.bam >test.bcf

From PL value, P(D|TT)=10^(-2.7), P(D|TA)=0. So TA allele should be the most possible allele.

But in the processed vcf using this command:

bcftools view -c test.bcf > test_count.vcf

This SNP is:


There is no alternative base information, but the AF1 is about 0.5, which means this is an heterozygous allele. And, if I use the following command, there is no SNP called at all.

bcftools view -vcg test.bcf > test_count.vcf

Why the alternative base is not A??? Is there any option I can choose to call this SNP?

BTW, I am using samtools version 0.1.10-6.

Cheers,

Jun

Still unsolved... waiting for an answer...

[lh3]

you need to apply the prior

[nkwuji]

Originally Posted by lh3

you need to apply the prior

Hi, Lh3,

Thx for your reply. Could you elaborate it? I find an option in bcttools view

-P STR, type of prior: full, cond2, flat

Do you mean this option? What are the differences between them?

Cheers,
Jun

[LisaM]

Originally Posted by nkwuji

Hi, Lh3,

Thx for your reply. Could you elaborate it? I find an option in bcttools view

-P STR, type of prior: full, cond2, flat

Do you mean this option? What are the differences between them?

Cheers,
Jun

Hi, could someone please comment on this please? I'm also looking into the differences between the "prior" commands, and I haven't been able to find anything on the internet. Any information would be helpful.

Thanks!
Lisa

[whodo_voodoo]

Hi,

Another request here for some more information with regards the 'prior' and snp calling options of bcftools. I've been trying to learn how to change the parameters of the snp calling but can't find much information on how to change the parameters in bcftools.

I've also got a second related question about bcftools, the mpileup page says

BCFtools does not properly handle multi-allelic variants. It only takes the strongest non-reference allele.

Does anybody know how much of an issue this is? I'm expecting a number of my SNPs to be heterozygous with no allele match to the reference, will this issue prevent me seeing them?

Thanks in advance for any answers,

Craig

[dan]

Originally Posted by lh3

you need to apply the prior

-- (minus minus)