Use of LazyRecord with queried IndexedReader (VCF), is it possible?

[davidmasp]

Disclaimer: I am relatively new to rust/noodles so I am might be missing super obvious from the documentation.

I am trying to use the LazyRecord (vcf) in a queried IndexedReader without much luck.

So far I can tell that the reader has a read_lazy_record which would give me a lazy record, but in my previous implementation I was using the method query and as far as I can tell that iterates over noodles::vcf::Record right?

Is there any way to iterate over lazy records instead? If not, is there any way to do something similar? I am thinking I should be able to "seek" the file but I am a bit lost with that part of the documentation.

Thanks a lot for the help! 🙌🏻

[zaeleus]

There's currently no way to use the VCF query iterator with lazy records. You can, however, copy and modify the existing implementation to work with them, e.g.,

vcf_232.rs

// cargo add [email protected] --features bgzf,csi,core,vcf

use std::{
    env,
    io::{self, Read, Seek},
    str,
};

use noodles::{
    bgzf,
    core::{region::Interval, Position, Region},
    csi::{self, binning_index::index::reference_sequence::bin::Chunk, BinningIndex},
    vcf,
};

fn main() -> Result<(), Box<dyn std::error::Error>> {
    let mut args = env::args().skip(1);
    let src = args.next().expect("missing src");
    let region = args.next().expect("missing region").parse()?;

    let mut reader = vcf::indexed_reader::Builder::default().build_from_path(src)?;
    let _header = reader.read_header()?;

    let query = reader_query_lazy_records(&mut reader, &region)?;

    for result in query {
        let record = result?;
        let end = record_end(&record)?;
        println!("{} {}-{}", record.chromosome(), record.position(), end);
    }

    Ok(())
}

fn reader_query_lazy_records<'r, R>(
    reader: &'r mut vcf::IndexedReader<R>,
    region: &Region,
) -> io::Result<QueryLazyRecords<'r, R>>
where
    R: Read + Seek,
{
    let index = reader.index();

    let (reference_sequence_id, reference_sequence_name) = resolve_region(index, region)?;
    let chunks = index.query(reference_sequence_id, region.interval())?;

    Ok(QueryLazyRecords::new(
        reader.get_mut(),
        chunks,
        reference_sequence_name,
        region.interval(),
    ))
}

fn resolve_region(index: &dyn BinningIndex, region: &Region) -> io::Result<(usize, Vec<u8>)> {
    let header = index
        .header()
        .ok_or_else(|| io::Error::new(io::ErrorKind::InvalidInput, "missing tabix header"))?;

    let region_name = str::from_utf8(region.name())
        .map_err(|e| io::Error::new(io::ErrorKind::InvalidInput, e))?;

    let i = header
        .reference_sequence_names()
        .get_index_of(region_name)
        .ok_or_else(|| {
            io::Error::new(
                io::ErrorKind::InvalidInput,
                format!(
                    "region reference sequence does not exist in reference sequences: {region:?}"
                ),
            )
        })?;

    Ok((i, region.name().into()))
}

pub struct QueryLazyRecords<'r, R>
where
    R: Read + Seek,
{
    reader: vcf::Reader<csi::io::Query<'r, R>>,
    reference_sequence_name: Vec<u8>,
    interval: Interval,
    record: vcf::lazy::Record,
}

impl<'r, R> QueryLazyRecords<'r, R>
where
    R: Read + Seek,
{
    fn new(
        reader: &'r mut bgzf::Reader<R>,
        chunks: Vec<Chunk>,
        reference_sequence_name: Vec<u8>,
        interval: Interval,
    ) -> Self {
        Self {
            reader: vcf::Reader::new(csi::io::Query::new(reader, chunks)),
            reference_sequence_name,
            interval,
            record: vcf::lazy::Record::default(),
        }
    }

    fn next_record(&mut self) -> io::Result<Option<vcf::lazy::Record>> {
        self.reader
            .read_lazy_record(&mut self.record)
            .map(|n| match n {
                0 => None,
                _ => Some(self.record.clone()),
            })
    }
}

impl<'r, R> Iterator for QueryLazyRecords<'r, R>
where
    R: Read + Seek,
{
    type Item = io::Result<vcf::lazy::Record>;

    fn next(&mut self) -> Option<Self::Item> {
        loop {
            match self.next_record() {
                Ok(Some(record)) => {
                    match intersects(&record, &self.reference_sequence_name, self.interval) {
                        Ok(true) => return Some(Ok(record)),
                        Ok(false) => {}
                        Err(e) => return Some(Err(e)),
                    }
                }
                Ok(None) => return None,
                Err(e) => return Some(Err(e)),
            }
        }
    }
}

pub(crate) fn intersects(
    record: &vcf::lazy::Record,
    reference_sequence_name: &[u8],
    region_interval: Interval,
) -> io::Result<bool> {
    let name = record.chromosome().to_string();

    let start = parse_position(record.position())?;
    let end = record_end(record)?;
    let record_interval = Interval::from(start..=end);

    Ok(name.as_bytes() == reference_sequence_name && record_interval.intersects(region_interval))
}

fn record_end(record: &vcf::lazy::Record) -> io::Result<Position> {
    const END: &str = "END";

    let start = parse_position(record.position())?;

    if let Some(Some(value)) = info_get(&record.info(), END) {
        parse_position(value)
    } else {
        let len = record.reference_bases().len();
        assert!(len > 0);
        Ok(start
            .checked_add(len - 1)
            .expect("attempt to add with overflow"))
    }
}

fn info_get<'a>(info: &'a vcf::lazy::record::Info<'a>, key: &str) -> Option<Option<&'a str>> {
    info.iter().find(|(k, _)| *k == key).map(|(_, v)| v)
}

fn parse_position(s: &str) -> io::Result<Position> {
    s.parse::<usize>()
        .map_err(|e| io::Error::new(io::ErrorKind::InvalidData, e))
        .and_then(|n| {
            Position::try_from(n).map_err(|e| io::Error::new(io::ErrorKind::InvalidData, e))
        })
}

Of course, this isn't ideal but is a usable workaround.

I'll note that the variant formats are undergoing a rework, and in the future, the query iterator will emit lazy records by default.

[davidmasp]

That makes sense, thanks for the answer and the work!

[zaeleus]

I'll note that the variant formats are undergoing a rework, and in the future, the query iterator will emit lazy records by default.

The query iterator in noodles 0.68.0 / noodles-vcf 0.52.0 now returns (lazy) format records by default.

[davidmasp]

Awesome! thanks for the update 👍🏻