samtools/bcftools strange (error ?) lack of call

Question

[mbourgey]

Hi all,

I got a strange lack of call when I used samtolls mpileup + bcftools to call my SNP (version 0.19 or 1.0 give me the same outputs)

Here as what I observe when I use samtools+bcftools view:

samtools mpileup -q 1 -u -D -S -g -f hg1k_v37.fasta -r 3:44775917-44775917 alignment/sample1.sorted.dup.recal.bam alignment/sample2.sorted.dup.recal.bam alignment/sample3.sorted.dup.recal.bam | bcftools view

##ALT=<ID=X,Description="Represents allele(s) other than observed.">
##INFO=<ID=INDEL,Number=0,Type=Flag,Description="Indicates that the variant is an INDEL.">
##INFO=<ID=IDV,Number=1,Type=Integer,Description="Maximum number of reads supporting an indel">
##INFO=<ID=IMF,Number=1,Type=Float,Description="Maximum fraction of reads supporting an indel">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Raw read depth">
##INFO=<ID=VDB,Number=1,Type=Float,Description="Variant Distance Bias for filtering splice-site artefacts in RNA-seq data (bigger is better)",Version=3>
##INFO=<ID=RPB,Number=1,Type=Float,Description="Mann-Whitney U test of Read Position Bias (bigger is better)">
##INFO=<ID=MQB,Number=1,Type=Float,Description="Mann-Whitney U test of Mapping Quality Bias (bigger is better)">
##INFO=<ID=BQB,Number=1,Type=Float,Description="Mann-Whitney U test of Base Quality Bias (bigger is better)">
##INFO=<ID=MQSB,Number=1,Type=Float,Description="Mann-Whitney U test of Mapping Quality vs Strand Bias (bigger is better)">
##INFO=<ID=SGB,Number=1,Type=Float,Description="Segregation based metric.">
##INFO=<ID=MQ0F,Number=1,Type=Float,Description="Fraction of MQ0 reads (smaller is better)">
##INFO=<ID=I16,Number=16,Type=Float,Description="Auxiliary tag used for calling, see description of bcf_callret1_t in bam2bcf.h">
##INFO=<ID=QS,Number=R,Type=Float,Description="Auxiliary tag used for calling">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="List of Phred-scaled genotype likelihoods">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Number of high-quality bases">
##FORMAT=<ID=SP,Number=1,Type=Integer,Description="Phred-scaled strand bias P-value">
##bcftools_viewVersion=1.0+htslib-1.0
##bcftools_viewCommand=view
#CHROM    POS    ID    REF    ALT    QUAL    FILTER    INFO    FORMAT    sample1    sample2    sample3
3    44775917    .    A    C,<X>    0    .    DP=247;I16=207,1,36,0,6174,186100,1117,35331,12480,748800,2160,129600,3824,84932,643,14149;QS=2.65514,0.344865,0;VDB=0.986901;SGB=0.623532;RPB=0.19008;MQB=1;MQSB=1;BQB=2.58832e-10;MQ0F=0    PL:DP:SP    0,90,152,90,152,152:30:0    20,0,115,255,185,236:112:0    0,42,110,255,149,215:102:0

when I then use: bcftools call -m -v I didn't get any variant called whereas the PL field of the sample 2 (20,0,115,255,185,236) show an higher likelihood for the heterozygote call.

Any idea what's happens here?

Thanks in advance

Mathieu

[Devon Ryan]

Hello mbourgey!

It appears that your post has been cross-posted to another site: samtools-help list

This is typically not recommended as it runs the risk of annoying people in both communities.

[mbourgey]

yes it 's me

sorry about that !

I posted here then think perhaps it should be more appropriate to ask that on the mailing list.

Mathieu

[Cytosine]

I don't believe it is an issue of genotype probabilites. What value did you feed to the option -m?

[mbourgey]

Hi Cytosine,

Thanks for your reply.

I didn't feed any value for the -m option in bcftools call as it does not require anyone:

-m, --multiallelic-caller alternative model for multiallelic and rare-variant calling (conflicts with -c)

and let the default value (not set) in samtools mpileup for the -m option:

-m, --min-ireads INT minimum number gapped reads for indel candidates [1]

Mathieu

[Cytosine]

Strange... I have bcftools v0.1.19-44428cd and the option menu says this:

-m FLOAT alternative model for multiallelic and rare-variant calling, include if P(chi^2)>=FLOAT

Meaning filter according to the chi^2 statistic... I thought your output was empty since you did not provide a value to -m

On the other hand, your option does tell something interesting (conflicts with -c)

And you call bcftools view also with -v, which makes it imply -c as well:

-v output potential variant sites only (force -c)

Could be that these two are interfering in your case, since -v forces -c, which conflicts with your -m.

[mbourgey]

Just for your information: it's the commands and options I used with samtools 1.0

I will check for the -c vs. -m conflict

Thanks
Mathieu

[mbourgey]

In bcftools 1.0 the -v option does not imply/force -c.

But just to be sure I tried using bcftools call -c -v instead and the result is the same no variant called