Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

consensus --missing not compatible with gvcf blocks #2350

Closed
Fan-iX opened this issue Jan 7, 2025 · 1 comment
Closed

consensus --missing not compatible with gvcf blocks #2350

Fan-iX opened this issue Jan 7, 2025 · 1 comment

Comments

@Fan-iX
Copy link

Fan-iX commented Jan 7, 2025

bcftools version: 1.21
Doing consensus with flag --missing will throw error on missing (./.) gvcf blocks. Here is a minimal reproducible example:

test.vcf.gz:

##fileformat=VCFv4.2
##contig=<ID=1,length=100>
##ALT=<ID=*,Description="Represents allele(s) other than observed.">
##INFO=<ID=END,Number=1,Type=Integer,Description="End position of the variant described in this record">
##INFO=<ID=MIN_DP,Number=1,Type=Integer,Description="Minimum per-sample depth in this gVCF block">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  1
1       3       .       A       .       .       .       END=5;MIN_DP=1;AN=0     GT      ./.
Inline commands to reproduce, for your convenience
printf '##fileformat=VCFv4.2\n##contig=<ID=1,length=100>\n##ALT=<ID=*,Description="Represents allele(s) other than observed.">\n##INFO=<ID=END,Number=1,Type=Integer,Description="End position of the variant described in this record">\n##INFO=<ID=MIN_DP,Number=1,Type=Integer,Description="Minimum per-sample depth in this gVCF block">\n##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">\n##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">\n#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t1\n1\t3\t.\tA\t.\t.\t.\tEND=5;MIN_DP=1;AN=0\tGT\t./.' | bcftools view -Oz -o test.vcf.gz --write-index

Running command

printf ">1\nAAAAAA" | bcftools consensus --missing N test.vcf.gz

will throws error:

The fasta sequence does not match the REF allele at 1:3:
   REF .vcf: [A]
   ALT .vcf: [N]
   REF .fa : [AAA]A

While the expected output is

>1
AANNNA

Or is there any other way to preserve missing sites when performing consensus on gvcf files?

@pd3 pd3 closed this as completed in 724713f Jan 25, 2025
@pd3
Copy link
Member

pd3 commented Jan 25, 2025

Thank you for the issue, this is now possible. Note, however, that the change will take effect only after samtools/htslib#1879 is merged

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants