Hello,

I've been running into an error lately when I try to run bcftools stats but that I believe is probably an error produced in bcftools call.

I have a set of .bcf files generated with samtools mpileup. Original bams were mapped with BWA-mem, treated with AddOrReplaceReadGroups, MarkDuplicates and RealignerTargetCreator.

This is my call for the files; I run them inside a lop

samtools mpileup -Q 20 -q 20 -go %s_raw.bcf -f %s %s

where the second %s is a path to the reference and the third the name of the input bam.

Then I run bcftools call, again inside a loop:

cat chromosome.list | xargs -I {} -n 1 -P 20 sh -c 'bcftools call -cO z -r {} --format-fields GQ -o %s/%s.{}.vcf.gz %s'

where the call would look something like this:

cat chromosome.list | xargs -I {} -n 1 -P 20 sh -c 'bcftools call -cO z -r {} --format-fields GQ -o LER_ALM_1580_raw/LER_ALM_1580_raw.{}.vcf.gz LER_ALM_1580_raw.bcf'

then I concatenate all the chromosome files using bcftools concat.

Finally, I try to run bcftools stats, like so (example from one sample):

bcftools stats -F /scratch/3/msf/Reference/LER_ALM1580/LER_ALM1580.gatk.iteration4.consensus.FINAL.fa -s - LER_ALM_1580_raw.vcf.gz > LER_ALM_1580_raw_v2.vcf.gz.stats

And get this error:

[E::bcf_calc_ac] Incorrect allele ("1") in LER_ALM_1580 at 1:461

If I look it up with bcftools view:

bcftools view -r 1:461 LER_ALM_1580_raw.vcf.gz

This is what I see:

CHROM POS ID REF ALT QUAL FILTER INFO FORMAT LER_ALM_1580

1 461 . C . 3.53256 . DP=1;SGB=-0.379885;MQ0F=0;AF1=1;AC1=2;DP4=0,0,0,1;MQ=60;FQ=-29.9956 GT:PL:GQ 0/1:29:5

Which means that there is a discordance between the GT field and REF/ALT information.

I have gone back and tried different things with only one sample and one chromosome. I've realized that if I repeat the analysis without using format-fields, the output vcf file is different (REF/ALT and GT are concordant!):

bcftools call -cO z -r 1 -o LER_ALM_1580_raw_chr1.vcf.gz LER_ALM_1580_raw.bcf
bcftools index LER_ALM_1580_raw_chr1.vcf.gz
bcftools stats -F /scratch/3/msf/SpeciesForReference/LER_ALM1580/LER_ALM1580.gatk.iteration4.consensus.FINAL.fa -s - LER_ALM_1580_raw_chr1.vcf.gz > LER_ALM_1580_raw_chr1.vcf.gz.stats
bcftools view -r 1:461 LER_ALM_1580_raw_chr1.vcf.gz

CHROM POS ID REF ALT QUAL FILTER INFO FORMAT LER_ALM_1580 1 461 . C . 3.53256 . DP=1;SGB=-0.379885;MQ0F=0;AF1=1;AC1=2;DP4=0,0,0,1;MQ=60;FQ=-29.9956 GT:PL 0/0:29

But if I use --format-fields, then the same as before occurs.

Any idea why this happens? Any advice on how to proceed?

Thank you