deserializing records from htsget stream

[cmdoret]

Dear developers,

When streaming chunks using noodles-htsget, is there a way to lazily parse CRAM records in memory?

In the noodles-htsget example, I've tried to replace:

noodles/noodles-htsget/examples/htsget_download_reads.rs

Lines 32 to 34 in 779500a

	while let Some(chunk) = chunks.try_next().await? {
	stdout.write_all(&chunk).await?;
	}

with things like:

let first_chunk = &chunks
    .try_next()
    .await
    .expect("can't load chunk")
    .unwrap()[..];

let mut reader = Reader::new(first_chunk);
let header = reader
    .read_header()
    .unwrap();

while let Some(chunk) = chunks.try_next().await.expect("can't load chunk") {
  let mut reader = Reader::new(&chunk[..]);
  for result in reader.records(&header) {
    let record = result.expect(&format!("{:?}", chunk));
      println!("{:?}", record);
    }
}

But when trying to instantiate records, I keep getting:

Custom { kind: InvalidData, error: TryFromInt
Error(()) }

Would you have any advice to parsing these records without saving the actual file?

[zaeleus]

I would recommend converting the chunks stream to an async reader using tokio_util::io::StreamReader. This way, you can use cram::r#async::io::Reader as you normally would.

Since you're working with CRAM, note that it commonly requires a reference sequence. Also note that htsget responses may include records outside the requested region(s), which you must filter yourself.

Here's a full example showing all of this:

main.rs

// cargo add [email protected] --no-default-features --features std
// cargo add [email protected] --features async,core,cram,fasta,htsget,sam
// cargo add [email protected] --features io-std,macros,rt-multi-thread
// cargo add [email protected] --features io

use std::env;

use futures::{StreamExt, TryStreamExt};
use noodles::{
    core::{region::Interval, Region},
    cram, fasta,
    htsget::{self, reads::Format},
    sam,
};
use tokio::io;
use tokio_util::io::StreamReader;

#[tokio::main]
async fn main() -> Result<(), Box<dyn std::error::Error>> {
    let mut args = env::args().skip(1);

    let base_url = args.next().expect("missing base URL").parse()?;
    let id = args.next().expect("missing ID");
    let fasta_src = args.next().expect("missing FASTA src");
    let region: Region = args.next().expect("missing region").parse()?;

    let reference_sequence_repository = fasta::io::indexed_reader::Builder::default()
        .build_from_path(fasta_src)
        .map(fasta::repository::adapters::IndexedReader::new)
        .map(fasta::Repository::new)?;

    let client = htsget::Client::new(base_url);

    let request = client
        .reads(id)
        .set_format(Format::Cram)
        .add_region(region.clone());

    let reads = request.send().await?;

    let chunks = reads
        .chunks()
        .map(|result| result.map_err(|e| io::Error::new(io::ErrorKind::InvalidData, e)));
    let chunks_reader = StreamReader::new(chunks);

    let mut reader = cram::r#async::io::reader::Builder::default()
        .set_reference_sequence_repository(reference_sequence_repository)
        .build_from_reader(chunks_reader);

    let header = reader.read_header().await?;

    let mut records = reader.records(&header);

    let mut writer = sam::r#async::io::Writer::new(io::stdout());
    writer.write_header(&header).await?;

    let reference_sequence_id = header
        .reference_sequences()
        .get_index_of(region.name())
        .expect("missing reference sequence");

    let region_interval = region.interval();

    while let Some(record) = records.try_next().await? {
        if !intersects(&record, reference_sequence_id, region_interval) {
            continue;
        }

        writer.write_alignment_record(&header, &record).await?;
    }

    Ok(())
}

fn intersects(record: &cram::Record, reference_sequence_id: usize, interval: Interval) -> bool {
    if let (Some(id), Some(start), Some(end)) = (
        record.reference_sequence_id(),
        record.alignment_start(),
        record.alignment_end(),
    ) {
        id == reference_sequence_id && interval.intersects((start..=end).into())
    } else {
        false
    }
}

[cmdoret]

Thank you for the comprehensive answer, this is really awesome!