search

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
640 stars 241 forks source link

`bcftools norm --multiallelics -indels` does not normalize indels. #2084

Closed ghuls closed 5 months ago

ghuls commented 6 months ago

bcftools norm --multiallelics -indels and bcftools norm --multiallelics +indels do not normalize indels even if all multiallelic sites are indels. with -snp and +snp it works for snps (altough if you have a snp and multiallelic indels, both snp and each indel gets a new entry like they would with -both).

❯  bcftools version
bcftools 1.19
Using htslib 1.19
Copyright (C) 2023 Genome Research Ltd.
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law.

# VCF file:
#  - first entry: multiallelic snps
#  - second entry: snp + multiallelic indels
#  - third entry: multiallelic indels
❯   bcftools view test_norm.vcf | bcftools view -H 
chr1    29291   .       C       T,G,A   10.85   PASS    F;AN=2;AC=4,3,4 GT:GQ:DP:AF     ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.     ./.:.:.:.        ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.      ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.     ./.:.:.:.        ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.      ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.     ./.:.:.:.        ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.      ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.     ./.:.:.:.        ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.      ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.     ./.:.:.:.        ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./3:.:.:.       ./.:.:.:.       ./1:.:.:.       ./.:.:.:.       ./.:.:.:.      ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./3:.:.:.       ./.:.:.:.       0/1:10:6:0.3333 ./.:.:.:.       ./1:.:.:.     ./1:.:.:.        ./2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./3:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./3:.:.:.
chr2    29292   .       T       C,TCCCTCTCCTTTCTCCTCTCTAGCC,TCTCTTTCTCACTGTCTCTCTAGCC,TCCCTCTCCTTTCTCCTCTCTAGC,TCCATCTGTATCCTCTCTAAGC,TCCCTCTCCTTTCTCCTCAGCC,TCCCTCTCCCTTTCTCCTCTCTAGCC,TCCTCTCCTTTCTCCTCTACCGC,TCCCTCTCCTTTCTCTCTCTAGCC,TCCCTCTCCTTTCTCCTCTAGCC,TCCCTCTCCTTTTCCTCCCCAGCC,TCCCTCTCCTTCTCCTCTCTAGCC,TCCCTCTCCCTTCTCCTCTCTCAC    14.66   PASS    F;AN=87;AC=9,57,1,3,2,4,2,1,1,1,1,1,1   GT:GQ:DP:AF     ./.:.:.:.       2/2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./2:.:.:.      0/1:6:5:0.4     ./.:.:.:.       ./3:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./4:.:.:.       ./.:.:.:.       ./.:.:.:.       5/5:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.     ./.:.:.:.        ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       2/2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.      ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./6:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./2:.:.:.       ./2:.:.:.       2/2:.:.:.     ./.:.:.:.        ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       2/2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       2/2:.:.:.       ./2:.:.:.       2/2:.:.:.       ./.:.:.:.       2/2:.:.:.       ./.:.:.:.       ./.:.:.:.       2/2:.:.:.      7/7:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       2/2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./4:.:.:.       2/2:.:.:.       ./.:.:.:.       ./4:.:.:.       ./8:.:.:.     ./.:.:.:.        ./.:.:.:.       ./2:.:.:.       ./.:.:.:.       2/2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.      ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./2:.:.:.       1/1:2:7:0.5714  ./2:.:.:.       ./2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./2:.:.:.       2/2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./2:.:.:.       ./.:.:.:.     ./.:.:.:.        ./2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./2:.:.:.       ./2:.:.:.       ./.:.:.:.       ./9:.:.:.       ./6:.:.:.       ./2:.:.:.       ./.:.:.:.       ./2:.:.:.       ./.:.:.:.       2/2:.:.:.       ./.:.:.:.       2/2:.:.:.      ./2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./2:.:.:.       ./.:.:.:.       ./2:.:.:.       1/1:5:3:1       ./10:.:.:.      ./.:.:.:.       0/1:6:6:0.5     ./2:.:.:.     ./.:.:.:.        ./.:.:.:.       ./11:.:.:.      ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       0/1:3:7:0.4286  ./.:.:.:.       ./12:.:.:.      ./.:.:.:.       ./.:.:.:.       ./2:.:.:.      ./2:.:.:.       ./.:.:.:.       ./2:.:.:.       ./2:.:.:.       ./2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./6:.:.:.       ./.:.:.:.       ./.:.:.:.       ./2:.:.:.       ./.:.:.:.       ./13:.:.:.      ./.:.:.:.       1/1:5:6:0.8333./6:.:.:.        ./2:.:.:.
chr3    29292   .       T       CGTA,TCCCTCTCCTTTCTCCTCTCTAGCC,TCTCTTTCTCACTGTCTCTCTAGCC,TCCCTCTCCTTTCTCCTCTCTAGC,TCCATCTGTATCCTCTCTAAGC,TCCCTCTCCTTTCTCCTCAGCC,TCCCTCTCCCTTTCTCCTCTCTAGCC,TCCTCTCCTTTCTCCTCTACCGC,TCCCTCTCCTTTCTCTCTCTAGCC,TCCCTCTCCTTTCTCCTCTAGCC,TCCCTCTCCTTTTCCTCCCCAGCC,TCCCTCTCCTTCTCCTCTCTAGCC,TCCCTCTCCCTTCTCCTCTCTCAC 14.66   PASS    F;AN=87;AC=9,57,1,3,2,4,2,1,1,1,1,1,1   GT:GQ:DP:AF     ./.:.:.:.       2/2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./2:.:.:.      0/1:6:5:0.4     ./.:.:.:.       ./3:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./4:.:.:.       ./.:.:.:.       ./.:.:.:.       5/5:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.     ./.:.:.:.        ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       2/2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.      ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./6:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./2:.:.:.       ./2:.:.:.       2/2:.:.:.     ./.:.:.:.        ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       2/2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       2/2:.:.:.       ./2:.:.:.       2/2:.:.:.       ./.:.:.:.       2/2:.:.:.       ./.:.:.:.       ./.:.:.:.       2/2:.:.:.      7/7:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       2/2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./4:.:.:.       2/2:.:.:.       ./.:.:.:.       ./4:.:.:.       ./8:.:.:.     ./.:.:.:.        ./.:.:.:.       ./2:.:.:.       ./.:.:.:.       2/2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.      ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./2:.:.:.       1/1:2:7:0.5714  ./2:.:.:.       ./2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./2:.:.:.       2/2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./2:.:.:.       ./.:.:.:.     ./.:.:.:.        ./2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./2:.:.:.       ./2:.:.:.       ./.:.:.:.       ./9:.:.:.       ./6:.:.:.       ./2:.:.:.       ./.:.:.:.       ./2:.:.:.       ./.:.:.:.       2/2:.:.:.       ./.:.:.:.       2/2:.:.:.      ./2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./2:.:.:.       ./.:.:.:.       ./2:.:.:.       1/1:5:3:1       ./10:.:.:.      ./.:.:.:.       0/1:6:6:0.5     ./2:.:.:.     ./.:.:.:.        ./.:.:.:.       ./11:.:.:.      ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       ./2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./.:.:.:.       0/1:3:7:0.4286  ./.:.:.:.       ./12:.:.:.      ./.:.:.:.       ./.:.:.:.       ./2:.:.:.      ./2:.:.:.       ./.:.:.:.       ./2:.:.:.       ./2:.:.:.       ./2:.:.:.       ./.:.:.:.       ./.:.:.:.       ./6:.:.:.       ./.:.:.:.       ./.:.:.:.       ./2:.:.:.       ./.:.:.:.       ./13:.:.:.      ./.:.:.:.       1/1:5:6:0.8333./6:.:.:.        ./2:.:.:.

# Number of current entries.
❯  bcftools view test_norm.vcf | bcftools view -H | wc -l
3

# Number indel entries.
❯  bcftools view test_norm.vcf | bcftools view --type indels | bcftools view -H | wc -l
2

# Split all multiallelic sites, and keep only indels.
❯  bcftools view test_norm.vcf | bcftools norm --multiallelics -both | bcftools view --type indels | bcftools view -H | wc -l
Lines   total/split/joined/realigned/skipped:   3/3/0/0/0
24

# Split all multiallelic sites, and keep only indels and recombine multiallelic indels
#    ==> not working, multiallelic indels are not merged to one entry.
❯  bcftools view test_norm.vcf | bcftools norm --multiallelics -both | bcftools view --type indels | bcftools norm --multiallelics +indels | bcftools view -H | wc -l
Lines   total/split/joined/realigned/skipped:   3/3/0/0/0
Lines   total/split/joined/realigned/skipped:   24/0/0/0/0
24

# Split all multiallelic sites, and keep only indels and recombine all multiallelic positions
#   ==> working, but it is actually only combining multiallelic indels in this case (as we filtered for indels before).
❯  bcftools view test_norm.vcf | bcftools norm --multiallelics -both | bcftools view --type indels | bcftools norm --multiallelics +both | bcftools view -H | wc -l
Lines   total/split/joined/realigned/skipped:   3/3/0/0/0
Lines   total/split/joined/realigned/skipped:   24/0/2/0/0
2

# Split all multiallelic indels, and keep only indels.
#   ==> not working, should at least have splitted entry 3 in multiple
❯  bcftools view test_norm.vcf | bcftools norm --multiallelics -indels | bcftools view --type indels | bcftools view -H | wc -l
Lines   total/split/joined/realigned/skipped:   3/0/0/0/0
2

# Split all multiallic positions, and keep only indels, combining all multiallic positions again (indels only due to filtering)
# and try tro split them on indels only.
#  ==> not working.
❯  bcftools view test_norm.vcf | bcftools norm --multiallelics -both | bcftools view --type indels | bcftools norm --multiallelics +both | bcftools norm --multiallelics -indels | bcftools view -H | wc -l
Lines   total/split/joined/realigned/skipped:   3/3/0/0/0
Lines   total/split/joined/realigned/skipped:   24/0/2/0/0
Lines   total/split/joined/realigned/skipped:   2/0/0/0/0
2

# Split all multiallelic sites, and keep only snps.
#   ==> works
❯   bcftools view test_norm.vcf | bcftools norm --multiallelics -both | bcftools view --type snps | bcftools view -H | wc -l
Lines   total/split/joined/realigned/skipped:   3/3/0/0/0
4

# Split all multiallelic SNP sites, and keep only snps.
#   ==> works
❯   bcftools view test_norm.vcf | bcftools norm --multiallelics -snps | bcftools view --type snps | bcftools view -H | wc -l
Lines   total/split/joined/realigned/skipped:   3/2/0/0/0
4

# Split all multiallelic SNP sites, and keep only snps and recombined multiallelic snps.
#   ==> works
❯  bcftools view test_norm.vcf | bcftools norm --multiallelics -snps | bcftools view --type snps | bcftools norm --multiallelics +snps | bcftools view -H | wc -l
Lines   total/split/joined/realigned/skipped:   3/2/0/0/0
Lines   total/split/joined/realigned/skipped:   4/0/1/0/0
2

# Split all multiallelic SNP sites, and count all entries.
#  ==> maybe semi broken. Entry 2 which had a snp and multiallelic indels, both snp and each indel gets a new entry.
❯  bcftools view test_norm.vcf | bcftools norm --multiallelics -snps | bcftools view -H | wc -l
Lines   total/split/joined/realigned/skipped:   3/2/0/0/0
17

test_norm.vcf.txt

pd3 commented 5 months ago

I can confirm that this is a problem, the code was written with SNPs vs anything else in mind. This needs to be improved / fixed.

pd3 commented 5 months ago

Thank you for the report. This has been improved, somewhat, with the new behavior and caveats described in the commit message of https://github.com/samtools/bcftools/commit/7a4f8010662003962e49280018f59ebaa87beb3a