How do I get a Genotype from a Sample?

[afontenot]

I'm iterating over the Records of a file and passing &Samples to a function to do some processing. Example code:

fn do_something(samples: &vcf::variant::record_buf::Samples) {
    for sample in samples.values() {
        let gtval = sample.get("GT").unwrap().unwrap();
    }
}

This would appear to get me what I'm looking for - gtval has type Value, but actually getting the Genotype out of this is proving difficult. For instance, attempting a match expression like let Genotype(gt) = gtval results in an error: cannot match against a tuple struct which contains private fields.

In general all the methods that do anything with a Sample just give you the generic Value object, so I don't see how to improve on this.

The documentation for methods on Sample is essentially non-existent: https://docs.rs/noodles/latest/noodles/vcf/variant/record_buf/samples/struct.Sample.html

As an aside, I'll mention that I got forced into this problem by trying to convert my use of noodles to async. Actually this isn't true, I had followed the docs in switching out records() with record_bufs(), but in fact this seems not to be necessary because the stream API also lets you iterate over records().

[zaeleus]

Record Samples are represented as a data frame.

In a single sample variant file, for example, the genotype (GT) series will only have a single Value::Genotype. Use Samples::select to get the series and then Series::get at index 0 to get the first sample's value.

Genotype is an abstraction over the serialized genotype form. (See § 1.6.2.8 "Genotype fields: GT (String)" (2024-06-28)). It doesn't return the encoded value as a string but position-phasing pairs ((Option<usize>, Phasing)).

All together:

use std::io;

use noodles_vcf::{
    self as vcf,
    variant::record::samples::{keys::key, series::Value, Series},
};

fn main() -> io::Result<()> {
    const DATA: &[u8] = b"\
##fileformat=VCFv4.3
#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tsample0
sq0\t1\t.\tA\t.\t.\tPASS\t.\tGT\t0/0
sq0\t2\t.\tC\t.\t.\tPASS\t.\tGT\t0|1
";

    let mut reader = vcf::io::Reader::new(DATA);
    let header = reader.read_header()?;

    for result in reader.records() {
        let record = result?;

        let samples = record.samples();
        let series = samples.select(key::GENOTYPE).expect("missing GT series");

        let value = series
            .get(&header, 0)
            .expect("missing sample")
            .expect("missing GT value")?;

        let Value::Genotype(genotype) = value else {
            panic!("invalid GT value");
        };

        for result in genotype.iter() {
            let allele = result?;
            let (position, phasing) = allele;
            print!("({position:?}, {phasing:?}) ");
        }

        println!();
    }

    Ok(())
}

[afontenot]

Thanks for the answer!

• So the correct way of iterating over samples is to first select the column you want and then use that, rather than going sample by sample and using the column, or columns, as needed? Or is it also valid to do it this way:

    for result in reader.records() {
        let record = result?;

        for sample in record.samples().iter() {
            let value = sample.get(&header, key::GENOTYPE).unwrap().unwrap().unwrap();
            let Value::Genotype(genotype) = value else {
                panic!("invalid GT value");
            };

            for result in genotype.iter() {
                let allele = result?;
                let (position, phasing) = allele;
                print!("({position:?}, {phasing:?}) ");
            }

            println!();
        }
    }

• It looks like you are matching on Genotype exactly the same way I was attempting to do. The issue is that if you're iterating over record_bufs() instead of records(), you have to match with a different Genotype:

use noodles_vcf::variant::record_buf::samples::sample::Value::Genotype

and then you also need to import the trait to allow iterating over the resulting object:

use noodles_vcf::variant::record::samples::series::value::Genotype as GenotypeTrait;

Which is messy to say the least. I don't think I was using record_bufs() for any particularly good reason, so I'm going to return to using records() as shown in your example.

[zaeleus]

Or is it also valid to do it this way:

Row-oriented (in your example) or columnar (in my example) are both valid approaches. They're just different ways to view data over the samples table.

The issue is that if you're iterating over record_bufs() instead of records(), you have to match with a different Genotype:

Like vcf::Record, vcf::variant::RecordBuf fully implements vcf::variant::Record. Value buffers can be converted to the shared series value (noodles::vcf::variant::record::samples::series::Value).

use noodles::vcf::variant::{
    record::samples::series::Value, record_buf::samples::sample::value::Value as ValueBuf,
};

let value = ValueBuf::Genotype("0|0".parse()?);

let Value::Genotype(_genotype) = (&value).into() else {
    panic!("invalid GT value");
}