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

bcf::Writer::write_vcf_record writes wrong FORMAT/GT #135

Closed
malthesr opened this issue Nov 22, 2022 · 1 comment
Closed

bcf::Writer::write_vcf_record writes wrong FORMAT/GT #135

malthesr opened this issue Nov 22, 2022 · 1 comment
Assignees
Labels
bug Something isn't working vcf

Comments

@malthesr
Copy link
Contributor

malthesr commented Nov 22, 2022

Consider the following code:

use std::io;

use bcf::header::StringMaps;
use noodles_bcf as bcf;

use noodles_vcf as vcf;
use vcf::header::{
    format::Key,
    record::value::{
        map::{Contig, Format},
        Map,
    },
};

fn main() -> Result<(), Box<dyn std::error::Error>> {
    let header = vcf::Header::builder()
        .add_contig(Map::<Contig>::new("sq0".parse()?))
        .add_sample_name("sample")
        .add_format(Map::<Format>::from(Key::Genotype))
        .add_format(Map::<Format>::from(Key::ConditionalGenotypeQuality))
        .build();

    let record: vcf::Record = "sq0\t1\t.\tA\t.\t.\t.\t.\tGT:GQ\t0|0:13".parse()?;

    // Write record as BCF: FORMAT/GT field is missing
    let mut writer = bcf::Writer::new(io::stdout());
    writer.write_file_format()?;
    writer.write_header(&header)?;
    writer.write_vcf_record(&header, &StringMaps::from(&header), &record)?;

    // Write record as VCF: Works fine
    // let mut writer = vcf::Writer::new(io::stdout());
    // writer.write_header(&header)?;
    // writer.write_record(&record)?;

    Ok(())
}

Running cargo run | bcftools view prints:

##fileformat=VCFv4.3
##FILTER=<ID=PASS,Description="All filters passed">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Conditional genotype quality">
##contig=<ID=sq0>
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	sample
sq0	1	.	A	.	.	.	.	GT:GQ	:13

Note the missing 0|0 genotype. bcftools also gives a warning [E::bcf_format_gt] Unexpected type 7.

noodles itself panics on trying to read the resulting BCF. For instance, running cargo run > $bcf and then running the bcf_view example on $bcf panics:

thread 'main' panicked at 'not yet implemented: unhandled type: Some(String(3))', noodles-bcf/src/reader/record/genotypes.rs:432:15
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

Writing the record as VCF seems to work fine, so I'm guessing the problem lies somewhere in bcf::Writer::write_vcf_record?

@zaeleus zaeleus added bug Something isn't working vcf labels Nov 22, 2022
@zaeleus zaeleus self-assigned this Nov 22, 2022
@zaeleus
Copy link
Owner

zaeleus commented Nov 28, 2022

Thanks for reporting and providing an example!

The BCF writer didn't have a FORMAT GT field encoder and was incorrectly writing string values instead. There is now support for this in 4c760ff.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working vcf
Projects
None yet
Development

No branches or pull requests

2 participants