bcftools icon indicating copy to clipboard operation
bcftools copied to clipboard

bcftools merge does not work

Open sjellerstrand opened this issue 1 year ago • 3 comments

Hej! I have found a potential bug in the develomental version. I get an error message when trying to merge several vcf files with bcftools merge. I have tried the the latest source code commit b7b2a32 for the following code. However, version 1.14, 1.16 and 1.17 works fine.

Command:

merge_files=$(echo ./Rasolark_2021_CADDXX010000136.1_HETGAM_95694_hetsex.vcf.gz ./Rasolark_2021_CADDXX010000136.1_HETGAM_TT95866_hetsex.vcf.gz)

bcftools merge $merge_files | bgzip -c > output.vcf.gz

The vcf file is not complete. It stops abruptly and has the following message at the end of the file:

[W::vcf_parse_info] INFO 'C' is not defined in the header, assuming Type=String
[E::bcf_write] Broken VCF record, the number of columns at CADDXX010000136.1:99689 does not match the number of samples (0 vs 2)
[main_vcfview] Error: cannot write to (null)

There is no INFO field called "C". What I do have is:

##INFO=<ID=AC,Number=A,Type=Integer,Description="Total number of alternate alleles in called genotypes">
##INFO=<ID=AF,Number=A,Type=Float,Description="Estimated allele frequency in the range (0,1]">
##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of samples with data">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">

The site in the individual vcf files look like this:

# Rasolark_2021_CADDXX010000136.1_HETGAM_95694_hetsex.vcf.gz
CADDXX010000136.1       99824   .       TA      T       0       .       AC=1;AF=1;AN=1;CM=0.091753;NS=18        GT      1

#Rasolark_2021_CADDXX010000136.1_HETGAM_TT95866_hetsex.vcf.gz
CADDXX010000136.1       99824   .       TA      T       0       .       AC=1;AF=1;AN=1;CM=0.091753;NS=18        GT      1

However, the problem does not seem to be that site in particular. The site changes depending on if I run more/other samples. For example if I merge 4 other samples the site 99824 is processed and the file instead ends at site 100730 with the message:

[E::bcf_write] Broken VCF record, the number of columns at CADDXX010000136.1:100730 does not match the number of samples (0 vs 4)
[main_vcfview] Error: cannot write to (null)

And if I merge all my 9 files I get the following message:

[E::vcf_parse_format] FORMAT column with no sample columns starting at CADDXX010000136.1:100636
Error: VCF parse error

Neither is there any obvious issue with this site in any of theindividual files:

#Rasolark_2021_CADDXX010000136.1_HETGAM_95694_hetsex.vcf.gz
CADDXX010000136.1       100636  .       C       G       0       .       AC=1;AF=0.25;AN=1;CM=0.092565;NS=18     GT      1
#Rasolark_2021_CADDXX010000136.1_HETGAM_TT95866_hetsex.vcf.gz
CADDXX010000136.1       100636  .       C       G       0       .       AC=1;AF=0.25;AN=1;CM=0.092565;NS=18     GT      1
#Rasolark_2021_CADDXX010000136.1_HETGAM_TT95867_hetsex.vcf.gz
CADDXX010000136.1       100636  .       C       G       0       .       AC=1;AF=0.25;AN=1;CM=0.092565;NS=18     GT      1
#Rasolark_2021_CADDXX010000136.1_HETGAM_TT95871_hetsex.vcf.gz
CADDXX010000136.1       100636  .       C       G       0       .       AC=1;AF=0.25;AN=1;CM=0.092565;NS=18     GT      1
#Rasolark_2021_CADDXX010000136.1_HETGAM_TT95879_hetsex.vcf.gz
CADDXX010000136.1       100636  .       C       G       0       .       AC=1;AF=0.25;AN=1;CM=0.092565;NS=18     GT      1
#Rasolark_2021_CADDXX010000136.1_HETGAM_TT95881_hetsex.vcf.gz
CADDXX010000136.1       100636  .       C       G       0       .       AC=1;AF=0.25;AN=1;CM=0.092565;NS=18     GT      1
#Rasolark_2021_CADDXX010000136.1_HETGAM_TT95884_hetsex.vcf.gz
CADDXX010000136.1       100636  .       C       G       0       .       AC=1;AF=0.25;AN=1;CM=0.092565;NS=18     GT      1
#Rasolark_2021_CADDXX010000136.1_HETGAM_TT95887_hetsex.vcf.gz
CADDXX010000136.1       100636  .       C       G       0       .       AC=1;AF=0.25;AN=1;CM=0.092565;NS=18     GT      1
#Rasolark_2021_CADDXX010000136.1_HETGAM_TT95888_hetsex.vcf.gz
CADDXX010000136.1       100636  .       C       G       0       .       AC=1;AF=0.25;AN=1;CM=0.092565;NS=18     GT      1

Both cases are fixed for the alternate allele, but other site like it are being processed without any issue.

Again, there is none of these problems in previous versions.

sjellerstrand avatar Jul 04 '23 12:07 sjellerstrand