bcftools loss record in targets when it comes with a G after pos.

[sssimonyang]

$ bcftools --version
bcftools 1.13
Using htslib 1.13

I have a sites.txt and a temp.vcf.

$ cat site.txt
1       10001   G       A
1       10002   A       A
1       10003   G       C
1       10004   A       C
$ cat temp.vcf
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##contig=<ID=1,length=20000000>
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO
1       10001   .       G       C       40      .       .
1       10002   .       T       A       50      .       .
1       10003   .       A       G       60      .       .
1       10004   .       A       T       70      .       .

Use bcftools view to get site in sites.txt.

$ bcftools view -T site.txt temp.vcf
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##contig=<ID=1,length=20000000>
##bcftools_viewVersion=1.13+htslib-1.13
##bcftools_viewCommand=view -T site.txt temp.vcf; Date=Mon Oct 18 20:03:23 2021
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO
1       10002   .       T       A       50      .       .
1       10004   .       A       T       70      .       .

The third columns of 10001 and 10003 is both G.

Only using the first two columns will help.

$ awk '{print $1"\t"$2}' site.txt > site1.txt
$ bcftools view -T site1.txt temp.vcf
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##contig=<ID=1,length=20000000>
##bcftools_viewVersion=1.13+htslib-1.13
##bcftools_viewCommand=view -T site1.txt temp.vcf; Date=Mon Oct 18 20:04:41 2021
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO
1       10001   .       G       C       40      .       .
1       10002   .       T       A       50      .       .
1       10003   .       A       G       60      .       .
1       10004   .       A       T       70      .       .

But it is unexpected that the unused columns have effect on result.

[pd3]

The program is trying to be smart and print only the matching REF and ALT positions. However, I just checked that this behavior is not documented and what is worse, it is not consistent: the matching is not attempted when the file is compressed and indexed with tabix.

[sssimonyang]

Please note that in above example, the 10002 and 10004 being queried out is not for their REF ALT being matched.
In site.txt, the record is 10004 A C and in vcf, it is 10004 A T.

So I think it is not due to matching REF and ALT.

[tobyaicher]

I have also just encountered an issue surrounding the -T flag for filtering and wanted to clarify.

1)I tried using the -T flag with two different site.txt files. One site.txt file had just CHROM and POS and the other had CHROM, POS, REF, and ALT

When I use this to filter a VCF file using -T, I notice that several variants that were present after selecting with CHROM and POS .txt were not present using the CHROM,POS, REF, and ALT .txt file. The sites missing after filtering with REF and ALT had the same REF and ALT in the vcf and site.txt file.

For all of these variants, the REF allele was a G.

I think there is a bug using a site.txt file that has a G as the REF base and that it causes bcftools to skip these sites.

2)Independent of the issue when the REF allele is a G, I am unsure whether or not bcftools view -T takes into account matching REF and ALT alleles or just the CHROM and POS sites. It seems like on the documentation suggests it is just the CHROM and POS sites? Is there a way to use -T and take into account the REF and ALT allele sites? This includes outside the context of a site.txt file (i.e. if I use -T using another vcf file to select the sites). If using -T only takes into account the CHROM and POS sites, is there a tool in bcftools that does take into REF and ALT? For example, for VCF file A only keep sites if the CHROM, POS, REF, and ALT are all identical in another VCF file B, but not if REF and/or ALT are different it will not keep this site.

Thank you let me know if I can clarify either of these points.

[pd3]

Thank you both for the report, apologies for taking so long to look at it. For now just want to acknowledge that the problem is indeed real and the cause has been traced and found to code in htslib.

[pd3]

This issue should resolve once samtools/htslib#1396 is merged.

[pd3]

This should be now resolved.