samtools / bcftools

This is the official development repository for BCFtools. See installation instructions and other documentation here http://samtools.github.io/bcftools/howtos/install.html
http://samtools.github.io/bcftools/
Other
652 stars 240 forks source link

bcftools call core dumped with -C alleles option #887

Closed pdroslva84 closed 5 years ago

pdroslva84 commented 5 years ago

Hi,

I'm trying to do genotype calling at known SNP positions, constraining the called alleles using the -C alleles option. My command line looks like this:

bcftools call -m -Oz -C alleles -T target_SNPs.txt.gz -o controls.vcf.gz controls.bcf

The controls.bcf file was generated using the same -T target_SNPs.txt.gz file with bcftools mpileup:

bcftools mpileup -Ob -o controls.bcf -f reference.fa sample1.sorted.bam sample2.sorted.bam -T target_SNPs.txt.gz -a AD,DP

However, with bcftools call I get the following error:

bcftools: mcall.c:1246: mcall_trim_numberR: Assertion `nret==nals*nsmpl' failed.
Aborted (core dumped)`

The problem seems to happen only if the target_SNPs.txt.gz file contains positions for which the bcf file does not contain the alt allele (because the samples genotyped happen to not have that allele). For example, for the following positions in the bcf:

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  sample1.sorted.bam       sample2.sorted.bam
chr1    212740  .       A       G,<*>   0       .       DP=73;I16=0,0,39,4,0,0,2743,192525,0,0,2580,154800,0,0,825,18621;QS=0,2,0;VDB=0.520868;SGB=-1.38232;MQSB=1;MQ0F=0       PL:DP:AD        255,72,0,255,72,255:24:0,24,0   255,57,0,255,57,255:19:0,19,0
chr1    320055  .       A       <*>     0       .       DP=101;I16=52,9,0,0,4116,300666,0,0,3660,219600,0,0,1281,29849,0,0;QS=2,0;MQSB=1;MQ0F=0 PL:DP:AD        0,87,255:29:29,0        0,96,255:32:32,0
chr1    486173  .       A       T,<*>   0       .       DP=13;I16=3,1,3,0,287,21853,246,20172,240,14400,180,10800,95,2275,75,1875;QS=1.25,0.75,0;VDB=0.074936;SGB=0.620439;RPB=0.810265;MQB=1.01283;MQSB=1;BQB=0.810265;MQ0F=0  PL:DP:AD        0,9,151,9,151,151:3:3,0,0      140,0,48,143,57,194:4:1,3,0
chr1    511277  .       A       G,<*>   0       .       DP=50;I16=0,0,25,4,0,0,1900,137374,0,0,1740,104400,0,0,672,16306;QS=0,2,0;VDB=0.0722735;SGB=-1.26186;MQSB=1;MQ0F=0      PL:DP:AD        255,30,0,255,30,255:10:0,10,0   255,57,0,255,57,255:19:0,19,0
chr1    602567  .       A       G,<*>   0       .       DP=9;I16=3,1,1,0,328,26896,41,1681,240,14400,60,3600,99,2451,18,324;QS=1.81448,0.18552,0;SGB=-0.516033;RPB=1;MQB=1;MQSB=1;BQB=1;MQ0F=0  PL:DP:AD        0,3,60,3,60,60:1:1,0,0  29,0,140,38,143,175:4:3,1,0
chr1    639707  .       T       A,<*>   0       .       DP=50;I16=0,0,23,8,0,0,1998,142356,0,0,1860,111600,0,0,612,13818;QS=0,2,0;VDB=0.563111;SGB=-1.37269;MQSB=1;MQ0F=0       PL:DP:AD        255,42,0,255,42,255:14:0,14,0   255,51,0,255,51,255:17:0,17,0

the following -T file works:

chr1    212740  A,G
chr1    486173  A,T
chr1    511277  A,G
chr1    602567  A,G
chr1    639707  T,A

but the following causes the crash (because of the 2nd position):

chr1    212740  A,G
chr1    320055  A,G
chr1    486173  A,T
chr1    511277  A,G
chr1    602567  A,G
chr1    639707  T,A

It works fine if not using the -C alleles option.

Any help would be appreciated. Thanks!

pd3 commented 5 years ago

Thank you for the bug report. The commit d88fa65 should fix it.

trhermes commented 2 years ago

Hello @pd3, are there plans to incorporate this commit into the main branch? We are having a similar issue as described above. Thanks!