noodles_sam/io/reader/
query.rs

1mod records;
2
3use std::io;
4
5use noodles_bgzf as bgzf;
6use noodles_core::region::Interval;
7use noodles_csi::{self as csi, binning_index::index::reference_sequence::bin::Chunk};
8
9use self::records::Records;
10use super::Reader;
11use crate::{Header, Record, alignment::Record as _};
12
13/// A reader over records of a SAM reader that intersects a given region.
14///
15/// This is created by calling [`Reader::query`].
16pub struct Query<'r, 'h: 'r, R> {
17    reader: Reader<csi::io::Query<'r, R>>,
18    header: &'h Header,
19    reference_sequence_id: usize,
20    interval: Interval,
21}
22
23impl<'r, 'h: 'r, R> Query<'r, 'h, R>
24where
25    R: bgzf::io::BufRead + bgzf::io::Seek,
26{
27    pub(super) fn new(
28        reader: &'r mut R,
29        chunks: Vec<Chunk>,
30        header: &'h Header,
31        reference_sequence_id: usize,
32        interval: Interval,
33    ) -> Self {
34        Self {
35            reader: Reader::new(csi::io::Query::new(reader, chunks)),
36            header,
37            reference_sequence_id,
38            interval,
39        }
40    }
41
42    /// Reads a record.
43    ///
44    /// # Example
45    ///
46    /// ```no_run
47    /// # use std::{fs::File, io};
48    /// use noodles_bgzf as bgzf;
49    /// use noodles_csi as csi;
50    /// use noodles_sam as sam;
51    ///
52    /// let mut reader = File::open("sample.sam.gz")
53    ///     .map(bgzf::io::Reader::new)
54    ///     .map(sam::io::Reader::new)?;
55    ///
56    /// let header = reader.read_header()?;
57    ///
58    /// let index = csi::fs::read("sample.sam.gz.csi")?;
59    /// let region = "sq0:8-13".parse()?;
60    /// let mut query = reader.query(&header, &index, &region)?;
61    ///
62    /// let mut record = sam::Record::default();
63    ///
64    /// while query.read_record(&mut record)? != 0 {
65    ///     // ...
66    /// }
67    /// # Ok::<(), Box<dyn std::error::Error>>(())
68    /// ```
69    pub fn read_record(&mut self, record: &mut Record) -> io::Result<usize> {
70        next_record(
71            &mut self.reader,
72            record,
73            self.header,
74            self.reference_sequence_id,
75            self.interval,
76        )
77    }
78
79    /// Returns an iterator over records.
80    pub fn records(self) -> Records<'r, 'h, R> {
81        Records::new(self)
82    }
83}
84
85pub(crate) fn intersects(
86    header: &Header,
87    record: &Record,
88    reference_sequence_id: usize,
89    region_interval: Interval,
90) -> io::Result<bool> {
91    let Some(id) = record.reference_sequence_id(header).transpose()? else {
92        return Ok(false);
93    };
94
95    if id != reference_sequence_id {
96        return Ok(false);
97    }
98
99    if interval_is_unbounded(region_interval) {
100        Ok(true)
101    } else {
102        match (
103            record.alignment_start().transpose()?,
104            record.alignment_end().transpose()?,
105        ) {
106            (Some(start), Some(end)) => {
107                let alignment_interval = (start..=end).into();
108                Ok(region_interval.intersects(alignment_interval))
109            }
110            _ => Ok(false),
111        }
112    }
113}
114
115fn interval_is_unbounded(interval: Interval) -> bool {
116    interval.start().is_none() && interval.end().is_none()
117}
118
119fn next_record<R>(
120    reader: &mut Reader<csi::io::Query<'_, R>>,
121    record: &mut Record,
122    header: &Header,
123    reference_sequence_id: usize,
124    interval: Interval,
125) -> io::Result<usize>
126where
127    R: bgzf::io::BufRead + bgzf::io::Seek,
128{
129    loop {
130        match reader.read_record(record)? {
131            0 => return Ok(0),
132            n => {
133                if intersects(header, record, reference_sequence_id, interval)? {
134                    return Ok(n);
135                }
136            }
137        }
138    }
139}
noodles_sam/io/reader/query.rs

noodles_sam/io/reader/
query.rs
