cli/commands/quantify: Filter unmapped records

Author: zaeleus

atlas-cli/src/commands/quantify.rs

@@ -1,4 +1,5 @@
 mod count;
+mod filter;
 mod match_intervals;
 mod specification;
 
@@ -18,7 +19,7 @@ use noodles::{bam, core::Position, gff::record::Strand, sam};
 use thiserror::Error;
 use tracing::info;
 
-use self::{count::count_single_records, specification::LibraryLayout};
+use self::{count::count_single_records, filter::Filter, specification::LibraryLayout};
 use crate::cli::quantify;
 
 type Features = HashMap<String, Vec<Feature>>;
@@ -72,11 +73,13 @@ pub fn quantify(args: quantify::Args) -> Result<(), QuantifyError> {
         "detected library layout"
     );
 
+    let filter = Filter;
+
     let mut reader = bam::io::reader::Builder.build_from_path(src)?;
     reader.read_header()?;
 
     let _ctx = match library_layout {
-        LibraryLayout::Single => count_single_records(&interval_trees, reader)?,
+        LibraryLayout::Single => count_single_records(&interval_trees, &filter, reader)?,
         LibraryLayout::Multiple => todo!(),
     };
 

atlas-cli/src/commands/quantify/count.rs

@@ -6,10 +6,10 @@ use std::{
 use atlas_core::collections::IntervalTree;
 use noodles::{bam, core::Position};
 
-use super::{match_intervals::MatchIntervals, Entry, IntervalTrees};
+use super::{match_intervals::MatchIntervals, Entry, Filter, IntervalTrees};
 
 #[expect(dead_code)]
-enum Event<'f> {
+pub(super) enum Event<'f> {
     Hit(&'f str),
     Miss,
     Ambiguous,
@@ -50,6 +50,7 @@ impl<'f> Context<'f> {
 
 pub(super) fn count_single_records<'f, R>(
     interval_trees: &IntervalTrees<'f>,
+    filter: &'f Filter,
     mut reader: bam::io::Reader<R>,
 ) -> io::Result<Context<'f>>
 where
@@ -59,7 +60,7 @@ where
     let mut record = bam::Record::default();
 
     while reader.read_record(&mut record)? != 0 {
-        let event = count_single_record(interval_trees, &record)?;
+        let event = count_single_record(interval_trees, filter, &record)?;
         ctx.add_event(event);
     }
 
@@ -68,8 +69,13 @@ where
 
 fn count_single_record<'f>(
     interval_trees: &IntervalTrees<'f>,
+    filter: &'f Filter,
     record: &bam::Record,
 ) -> io::Result<Event<'f>> {
+    if let Some(event) = filter.filter(record)? {
+        return Ok(event);
+    }
+
     let reference_sequence_id = record
         .reference_sequence_id()
         .transpose()?

atlas-cli/src/commands/quantify/filter.rs

@@ -0,0 +1,19 @@
+use std::io;
+
+use noodles::bam;
+
+use super::count::Event;
+
+pub(super) struct Filter;
+
+impl Filter {
+    pub(super) fn filter(&self, record: &bam::Record) -> io::Result<Option<Event<'_>>> {
+        let flags = record.flags();
+
+        if flags.is_unmapped() {
+            return Ok(Some(Event::Unmapped));
+        }
+
+        Ok(None)
+    }
+}