Handle CRAM references, empty thresholds, and zero-target extract inputs

modkit-core/Cargo.toml

@@ -29,6 +29,7 @@ indexmap = "2.2.6"
 indicatif = { version = "0.17.1", features = ["rayon"] }
 itertools = "0.12.1"
 lazy_static = "1.4"
+libc = "0.2"
 linear-map = "1.2.0"
 log = "0.4.0"
 log-once = "0.4.0"

modkit-core/src/command_utils.rs

@@ -13,7 +13,7 @@ use crate::mod_base_code::{DnaBase, ModCodeRepr};
 use crate::motifs::motif_bed::RegexMotif;
 use crate::position_filter::StrandedPositionFilter;
 use crate::threshold_mod_caller::MultipleThresholdModCaller;
-use crate::thresholds::calc_threshold_from_bam;
+use crate::thresholds::calc_threshold_from_bam_with_reference;
 use crate::util::{create_out_directory, Region};
 
 pub fn parse_per_mod_thresholds(
@@ -141,6 +141,44 @@ pub fn get_threshold_from_options(
     position_filter: Option<&StrandedPositionFilter<()>>,
     only_mapped: bool,
     suppress_progress: bool,
+) -> anyhow::Result<MultipleThresholdModCaller> {
+    get_threshold_from_options_with_reference(
+        in_bam,
+        None,
+        threads,
+        interval_size,
+        sample_frac,
+        num_reads,
+        no_filtering,
+        filter_percentile,
+        seed,
+        region,
+        per_mod_thresholds,
+        edge_filter,
+        collapse_method,
+        position_filter,
+        only_mapped,
+        suppress_progress,
+    )
+}
+
+pub(crate) fn get_threshold_from_options_with_reference(
+    in_bam: &PathBuf,
+    reference_fasta: Option<&PathBuf>,
+    threads: usize,
+    interval_size: u32,
+    sample_frac: Option<f64>,
+    num_reads: usize,
+    no_filtering: bool,
+    filter_percentile: f32,
+    seed: Option<u64>,
+    region: Option<&Region>,
+    per_mod_thresholds: Option<HashMap<ModCodeRepr, f32>>,
+    edge_filter: Option<&EdgeFilter>,
+    collapse_method: Option<&CollapseMethod>,
+    position_filter: Option<&StrandedPositionFilter<()>>,
+    only_mapped: bool,
+    suppress_progress: bool,
 ) -> anyhow::Result<MultipleThresholdModCaller> {
     if no_filtering {
         info!("not performing filtering");
@@ -157,8 +195,9 @@ pub fn get_threshold_from_options(
             (None, Some(num_reads))
         }
     };
-    let per_base_thresholds = calc_threshold_from_bam(
+    let per_base_thresholds = calc_threshold_from_bam_with_reference(
         in_bam,
+        reference_fasta,
         threads,
         interval_size,
         sample_frac,

modkit-core/src/entropy/mod.rs

@@ -25,7 +25,10 @@ use crate::read_ids_to_base_mod_probs::{PositionModCalls, ReadBaseModProfile};
 use crate::reads_sampler::sampling_schedule::ReferenceSequencesLookup;
 use crate::threshold_mod_caller::MultipleThresholdModCaller;
 use crate::thresholds::percentile_linear_interp;
-use crate::util::{record_is_not_primary, ReferenceRecord, Strand};
+use crate::util::{
+    record_is_not_primary, set_reference_for_cram_indexed_reader,
+    ReferenceRecord, Strand,
+};
 
 mod methylation_entropy;
 pub mod subcommand;
@@ -1486,11 +1489,13 @@ struct Message {
 
 fn process_bam_fp(
     bam_fp: &PathBuf,
+    reference_fasta: &PathBuf,
     fetch_definition: FetchDefinition,
     caller: Arc<MultipleThresholdModCaller>,
     io_threads: usize,
 ) -> anyhow::Result<Vec<Message>> {
     let mut reader = bam::IndexedReader::from_path(bam_fp)?;
+    set_reference_for_cram_indexed_reader(&mut reader, Some(reference_fasta))?;
     reader.set_threads(io_threads)?;
     reader.fetch(fetch_definition)?;
 
@@ -1583,9 +1588,11 @@ pub(super) fn process_entropy_window(
     io_threads: usize,
     caller: Arc<MultipleThresholdModCaller>,
     bam_fps: &[PathBuf],
+    reference_fasta: &PathBuf,
 ) -> anyhow::Result<EntropyCalculation> {
     let bam_fp = &bam_fps[0];
-    let reader = bam::IndexedReader::from_path(bam_fp)?;
+    let mut reader = bam::IndexedReader::from_path(bam_fp)?;
+    set_reference_for_cram_indexed_reader(&mut reader, Some(reference_fasta))?;
     let chrom_id = entropy_windows.chrom_id;
     drop(reader);
 
@@ -1594,6 +1601,7 @@ pub(super) fn process_entropy_window(
         .map(|fp| {
             process_bam_fp(
                 fp,
+                reference_fasta,
                 entropy_windows.get_fetch_definition(),
                 caller.clone(),
                 io_threads,

modkit-core/src/entropy/subcommand.rs

@@ -22,8 +22,8 @@ use crate::reads_sampler::sampling_schedule::{
 };
 use crate::threshold_mod_caller::MultipleThresholdModCaller;
 use crate::thresholds::{
-    calculate_threshold_with_fallback, get_modbase_probs_from_bam,
-    log_calculated_thresholds,
+    calculate_threshold_with_fallback,
+    get_modbase_probs_from_bam_with_reference, log_calculated_thresholds,
 };
 use crate::util::{
     format_errors_table, get_master_progress_bar, get_ticker, MutOpMax,
@@ -60,7 +60,12 @@ pub struct MethylationEntropy {
     window_size: usize,
     /// Do not perform any filtering, include all mod base calls in output.
     #[clap(help_heading = "Filtering Options")]
-    #[arg(group = "thresholds", long, default_value_t = false)]
+    #[arg(
+        group = "thresholds",
+        long,
+        conflicts_with = "mod_thresholds",
+        default_value_t = false
+    )]
     no_filtering: bool,
     /// Sample this many reads when estimating the filtering threshold. Reads
     /// will be sampled evenly across aligned genome. If a region is
@@ -194,42 +199,20 @@ impl MethylationEntropy {
             bail!("min-valid-coverage must be at least 1")
         }
         for bam_fp in self.in_bams.iter() {
-            IdxStats::check_any_mapped_reads(&bam_fp, None, None)
-                .with_context(|| {
-                    format!(
-                        "did not find any mapped reads in {bam_fp:?}, perform \
+            IdxStats::check_any_mapped_reads_with_reference(
+                &bam_fp,
+                Some(&self.reference_fasta),
+                None,
+                None,
+            )
+            .with_context(|| {
+                format!(
+                    "did not find any mapped reads in {bam_fp:?}, perform \
                          alignment first"
-                    )
-                })?;
+                )
+            })?;
         }
 
-        let mut writer: Box<dyn EntropyWriter> =
-            match (self.out_bed.as_ref(), self.regions_fp.is_some()) {
-                (Some(out_fp), false) => Box::new(
-                    WindowsWriter::new_file(out_fp, self.header, self.verbose)
-                        .context("failed to make writer to file")?,
-                ),
-                (Some(out_dir), true) => Box::new(
-                    RegionsWriter::new(
-                        out_dir,
-                        self.prefix.as_ref(),
-                        self.header,
-                        self.verbose,
-                    )
-                    .context(
-                        "failed to make regions writer, output must be a \
-                         directory",
-                    )?,
-                ),
-                (None, false) => Box::new(
-                    WindowsWriter::new_stdout(self.header, self.verbose)
-                        .context("failed to make writer to stdout")?,
-                ),
-                (None, true) => {
-                    bail!("must provide output directory with regions")
-                }
-            };
-
         let pool = rayon::ThreadPoolBuilder::new()
             .num_threads(self.threads)
             .build()?;
@@ -338,9 +321,43 @@ impl MethylationEntropy {
         let threshold_caller =
             self.get_threshold_caller(&pool).map(|c| Arc::new(c))?;
 
+        let mut writer: Box<dyn EntropyWriter> =
+            match (self.out_bed.as_ref(), self.regions_fp.is_some()) {
+                (Some(out_fp), false) => Box::new(
+                    WindowsWriter::new_file(
+                        out_fp,
+                        self.header,
+                        self.verbose,
+                        self.force,
+                    )
+                    .context("failed to make writer to file")?,
+                ),
+                (Some(out_dir), true) => Box::new(
+                    RegionsWriter::new(
+                        out_dir,
+                        self.prefix.as_ref(),
+                        self.header,
+                        self.verbose,
+                        self.force,
+                    )
+                    .context(
+                        "failed to make regions writer, output must be a \
+                         directory",
+                    )?,
+                ),
+                (None, false) => Box::new(
+                    WindowsWriter::new_stdout(self.header, self.verbose)
+                        .context("failed to make writer to stdout")?,
+                ),
+                (None, true) => {
+                    bail!("must provide output directory with regions")
+                }
+            };
+
         let (snd, rcv) = crossbeam::channel::bounded(10_000);
 
         let bam_fps = self.in_bams.clone();
+        let reference_fasta = self.reference_fasta.clone();
         let min_coverage = self.min_valid_coverage;
         let threads = self.threads;
         let io_threads = self.io_threads.unwrap_or(threads);
@@ -387,6 +404,7 @@ impl MethylationEntropy {
                                     io_threads,
                                     threshold_caller.clone(),
                                     &bam_fps,
+                                    &reference_fasta,
                                 )
                             })
                             .collect::<Vec<_>>();
@@ -455,6 +473,11 @@ impl MethylationEntropy {
         &self,
         pool: &rayon::ThreadPool,
     ) -> anyhow::Result<MultipleThresholdModCaller> {
+        if self.no_filtering {
+            info!("not performing filtering");
+            return Ok(MultipleThresholdModCaller::new_passthrough());
+        }
+
         let per_mod_thresholds = self
             .mod_thresholds
             .as_ref()
@@ -482,8 +505,9 @@ impl MethylationEntropy {
                     HashMap::<DnaBase, f32>::new();
                 for in_bam in self.in_bams.iter() {
                     let (per_base_thresholds, explicit_canonical_probs) =
-                        get_modbase_probs_from_bam(
+                        get_modbase_probs_from_bam_with_reference(
                             in_bam,
+                            Some(&self.reference_fasta),
                             self.threads,
                             1_000_000,
                             None,