FaiReader TabixReader

Cargo.toml

@@ -60,6 +60,7 @@ triple_accel = "0.4.0"
 chainfile = "0.3.0"
 omics = "0.2.0"
 filetime = "0.2.27"
+pandora_lib_smb = { git = "https://git.t0m4.fr/Thomas/pandora_lib_smb" }
 
 [profile.dev]
 opt-level = 0

src/callers/nanomonsv.rs

@@ -84,15 +84,11 @@ use crate::{
     annotation::{Annotation, Annotations, Caller, CallerCat, Sample},
     collection::vcf::Vcf,
     commands::{
-        bcftools::{BcftoolsConcat, BcftoolsKeepPass},
-        CapturedOutput, Command as JobCommand, LocalRunner, SbatchRunner, SlurmParams,
+        CapturedOutput, Command as JobCommand, LocalRunner, SbatchRunner, SlurmParams, bcftools::{BcftoolsConcat, BcftoolsKeepPass}
     },
     config::Config,
     helpers::{is_file_older, remove_dir_if_exists},
-    io::{
-        fasta::{open_indexed_fasta, reference_alternative_from_svlen},
-        vcf::read_vcf,
-    },
+    io::{fasta::{FaiReader, reference_alternative_from_svlen}, vcf::read_vcf},
     locker::SampleLock,
     pipes::{Initialize, InitializeSolo, ShouldRun, Version},
     run,
@@ -410,9 +406,9 @@ impl Version for NanomonSV {
     }
 }
 
-fn nanomonsv_var_norm(
+fn nanomonsv_var_norm<R: std::io::Read + std::io::Seek>(
     mut v: VcfVariant,
-    fa_reader: &mut noodles_fasta::io::IndexedReader<noodles_fasta::io::BufReader<std::fs::File>>,
+    fa_reader: &mut FaiReader<R>,
 ) -> VcfVariant {
     match v.svtype() {
         Some(SVType::DEL) => {
@@ -458,7 +454,8 @@ impl Variants for NanomonSV {
 
         let reference = self.config.reference.to_string();
         let path = PathBuf::from(reference);
-        let mut fa_reader = open_indexed_fasta(&path)?;
+        // let mut fa_reader = open_indexed_fasta_from_reader(&path)?;
+        let mut fa_reader = FaiReader::open(&path)?;
         let variants: Vec<VcfVariant> = read_vcf(&vcf_passed)
             .map_err(|e| anyhow::anyhow!("Failed to read NanomonSV VCF {vcf_passed}.\n{e}"))?
             .into_iter()
@@ -676,7 +673,7 @@ impl Variants for NanomonSVSolo {
 
         let reference = self.config.reference.to_string();
         let path = PathBuf::from(reference);
-        let mut fa_reader = open_indexed_fasta(&path)?;
+        let mut fa_reader = FaiReader::open(&path)?;
         let variants: Vec<VcfVariant> = read_vcf(&self.vcf_passed)
             .map_err(|e| anyhow::anyhow!("Failed to read NanomonSV VCF {}.\n{e}", self.vcf_passed))?
             .into_iter()

src/io/bed.rs

@@ -30,15 +30,14 @@ use rayon::prelude::*;
 
 use crate::{
     annotation::Annotation,
-    io::writers::{BgzTabixWriter, IndexFormat},
+    io::{readers::TabixReader, writers::{BgzTabixWriter, IndexFormat}},
     positions::{
-        extract_contig_indices, find_contig_indices, overlaps_par, GenomePosition, GenomeRange,
-        GetGenomeRange,
+        GenomePosition, GenomeRange, GetGenomeRange, extract_contig_indices, find_contig_indices, overlaps_par
     },
     variant::variant_collection::Variants,
 };
 
-use super::readers::{fetch_tabix_lines_with, get_reader};
+use super::readers::get_reader;
 
 /// One parsed row from a BED file (up to BED6).
 ///
@@ -413,9 +412,12 @@ pub fn convert_bgz_with_tabix(input: impl AsRef<Path>, force: bool) -> anyhow::R
 ///
 /// Returns an error if the index or file cannot be read, or if a fetched line
 /// fails to parse as a BED row.
-pub fn fetch_bed(bgz_path: &str, region: &GenomeRange) -> anyhow::Result<Vec<BedRow>> {
+pub fn fetch_bed_with_reader<R: std::io::Read + std::io::Seek>(
+    reader: &mut TabixReader<R>,
+    region: &GenomeRange,
+) -> anyhow::Result<Vec<BedRow>> {
     let mut rows = Vec::new();
-    fetch_tabix_lines_with(bgz_path, region, |line| {
+    reader.fetch_with(region, |line| {
         let row: BedRow = line
             .parse()
             .with_context(|| format!("failed to parse BED line: {line}"))?;

src/io/fasta.rs

@@ -6,6 +6,7 @@
 
 use std::collections::HashMap;
 use std::io::{BufRead, Write};
+use std::sync::Arc;
 use std::{
     fs::File,
     io::BufReader,
@@ -13,21 +14,76 @@ use std::{
 };
 
 use anyhow::Context;
-use noodles_fasta::io::{BufReader as NoodlesBufReader, IndexedReader};
+use noodles_fasta::{self as fasta, fai, io::indexed_reader::Builder, io::IndexedReader};
+use pandora_lib_smb::{SmbReader, SmbSession};
 
 use crate::{positions::GenomePosition, variant::vcf_variant::ReferenceAlternative};
 
-/// Open an indexed FASTA file (requires a `.fai` index alongside the path).
-///
-/// # Errors
-///
-/// Returns an error if the FASTA or its index cannot be opened.
-pub fn open_indexed_fasta(
-    reference: &Path,
-) -> anyhow::Result<IndexedReader<NoodlesBufReader<File>>> {
-    noodles_fasta::io::indexed_reader::Builder::default()
-        .build_from_path(reference)
-        .map_err(|e| anyhow::anyhow!("Failed to open indexed FASTA {}: {e}", reference.display()))
+
+pub struct FaiReader<R> {
+    inner: IndexedReader<std::io::BufReader<R>>,
+}
+
+impl<R> FaiReader<R>
+where
+    R: std::io::Read + std::io::Seek,
+{
+    /// Build from an already-open data reader plus pre-parsed `.fai` bytes.
+    fn from_reader(data: R, fai_bytes: &[u8]) -> anyhow::Result<Self> {
+        let index = fai::io::Reader::new(fai_bytes)
+            .read_index()
+            .context("parse .fai index")?;
+        let inner = Builder::default()
+            .set_index(index)
+            .build_from_reader(std::io::BufReader::new(data))
+            .map_err(|e| anyhow::anyhow!("Failed to build indexed FASTA: {e}"))?;
+        Ok(Self { inner })
+    }
+
+    /// e.g. delegate region queries — same shape as your existing usage.
+    pub fn query(&mut self, region: &noodles_core::Region) -> std::io::Result<fasta::Record> {
+        self.inner.query(region)
+    }
+}
+
+impl FaiReader<File> {
+    /// Open a local indexed FASTA from disk (replaces `open_indexed_fasta`).
+    pub fn open(reference: &Path) -> anyhow::Result<Self> {
+        let fai_path = reference.with_extension(
+            format!("{}.fai", reference.extension().and_then(|e| e.to_str()).unwrap_or("")),
+        );
+        // or just reference.to_string_lossy() + ".fai", matching your existing convention
+        let fai_bytes = std::fs::read(&fai_path)
+            .with_context(|| format!("read {}", fai_path.display()))?;
+        let data = File::open(reference)
+            .with_context(|| format!("open {}", reference.display()))?;
+        Self::from_reader(data, &fai_bytes)
+    }
+}
+
+impl FaiReader<SmbReader> {
+    /// Open an indexed FASTA hosted over SMB.
+    pub async fn open_smb(session: &Arc<SmbSession>, fasta_rel_path: &str) -> anyhow::Result<Self> {
+        let fai_rel_path = format!("{fasta_rel_path}.fai");
+        let fai_bytes = session
+            .download_file_bytes(&fai_rel_path)
+            .await
+            .context("download .fai over SMB")?;
+        let data = session
+            .open_smb_reader(fasta_rel_path)
+            .await
+            .context("open SMB reader for FASTA data")?;
+        Self::from_reader(data, &fai_bytes)
+    }
+}
+
+impl<R> std::ops::Deref for FaiReader<R> {
+    type Target = IndexedReader<std::io::BufReader<R>>;
+    fn deref(&self) -> &Self::Target { &self.inner }
+}
+
+impl<R> std::ops::DerefMut for FaiReader<R> {
+    fn deref_mut(&mut self) -> &mut Self::Target { &mut self.inner }
 }
 
 /// Fetch a fixed-width window of reference sequence centred on `pos0`.
@@ -40,8 +96,8 @@ pub fn open_indexed_fasta(
 ///
 /// Returns an error if the region is out of bounds or the sequence cannot be
 /// decoded as UTF-8.
-pub fn sequence_at(
-    fasta_reader: &mut IndexedReader<noodles_fasta::io::BufReader<File>>,
+pub fn sequence_at<R: std::io::Read + std::io::Seek>(
+    fasta_reader: &mut FaiReader<R>,
     contig: &str,
     pos0: usize,
     len: usize,
@@ -75,8 +131,8 @@ pub fn sequence_at(
 ///
 /// Returns an error if the region is out of bounds or the sequence cannot be
 /// decoded as UTF-8.
-pub fn sequence_range(
-    fasta_reader: &mut IndexedReader<noodles_fasta::io::BufReader<File>>,
+pub fn sequence_range<R: std::io::Read + std::io::Seek>(
+    fasta_reader: &mut FaiReader<R>,
     contig: &str,
     start0: usize,         // 0-based inclusive
     end0_inclusive: usize, // 0-based inclusive
@@ -106,8 +162,8 @@ pub fn sequence_range(
 ///
 /// This mirrors VCF's anchored allele representation, where indels include the
 /// base immediately before the changed sequence.
-pub fn reference_alternative_from_svlen(
-    fasta_reader: &mut IndexedReader<noodles_fasta::io::BufReader<File>>,
+pub fn reference_alternative_from_svlen<R: std::io::Read + std::io::Seek>(
+    fasta_reader: &mut FaiReader<R>,
     position: &GenomePosition,
     svlen: i32,
 ) -> anyhow::Result<ReferenceAlternative> {
@@ -124,8 +180,8 @@ pub fn reference_alternative_from_svlen(
     };
 
     sequence_range(fasta_reader, &position.contig(), start0, end0_inclusive)?
-    .parse()
-    .context("failed to parse FASTA sequence as ReferenceAlternative")
+        .parse()
+        .context("failed to parse FASTA sequence as ReferenceAlternative")
 }
 
 /// Fetch the VCF reference allele for an insertion from FASTA.
@@ -133,8 +189,8 @@ pub fn reference_alternative_from_svlen(
 /// Insertions are represented with an anchor base in VCF, so the reference
 /// allele is just the base at `position`. `svlen` is validated for consistency
 /// but cannot identify the inserted sequence by itself.
-pub fn insertion_reference_from_svlen(
-    fasta_reader: &mut IndexedReader<noodles_fasta::io::BufReader<File>>,
+pub fn insertion_reference_from_svlen<R: std::io::Read + std::io::Seek>(
+    fasta_reader: &mut FaiReader<R>,
     position: &GenomePosition,
     svlen: i32,
 ) -> anyhow::Result<ReferenceAlternative> {
@@ -150,8 +206,8 @@ pub fn insertion_reference_from_svlen(
 ///
 /// This is the correct helper for ordinary insertions where the inserted
 /// sequence comes from the caller (`SVINSSEQ`, sequence-resolved ALT, etc.).
-pub fn insertion_alternative_from_sequence(
-    fasta_reader: &mut IndexedReader<noodles_fasta::io::BufReader<File>>,
+pub fn insertion_alternative_from_sequence<R: std::io::Read + std::io::Seek>(
+    fasta_reader: &mut FaiReader<R>,
     position: &GenomePosition,
     inserted_sequence: &str,
 ) -> anyhow::Result<ReferenceAlternative> {
@@ -172,8 +228,8 @@ pub fn insertion_alternative_from_sequence(
 /// This is only valid for tandem-duplication-like insertions where the inserted
 /// bases are copied from the reference immediately after the anchor. For novel
 /// inserted sequence, use [`insertion_alternative_from_sequence`] instead.
-pub fn reference_derived_insertion_alternative_from_svlen(
-    fasta_reader: &mut IndexedReader<noodles_fasta::io::BufReader<File>>,
+pub fn reference_derived_insertion_alternative_from_svlen<R: std::io::Read + std::io::Seek>(
+    fasta_reader: &mut FaiReader<R>,
     position: &GenomePosition,
     svlen: i32,
 ) -> anyhow::Result<ReferenceAlternative> {
@@ -185,14 +241,9 @@ pub fn reference_derived_insertion_alternative_from_svlen(
         .checked_add(len)
         .context("SVLEN FASTA range end coordinate overflow")?;
 
-    sequence_range(
-        fasta_reader,
-        &position.contig(),
-        start0,
-        end0_inclusive,
-    )?
-    .parse()
-    .context("failed to parse FASTA sequence as ReferenceAlternative")
+    sequence_range(fasta_reader, &position.contig(), start0, end0_inclusive)?
+        .parse()
+        .context("failed to parse FASTA sequence as ReferenceAlternative")
 }
 
 /// A single-contig FASTA file produced by [`split_fasta`].

src/io/modkit.rs

@@ -35,8 +35,8 @@ use std::{
 };
 
 use crate::{
-    helpers::{diverging_rgb, TempFileGuard},
-    io::{bed::read_bed, readers::fetch_tabix_lines_with},
+    helpers::{TempFileGuard, diverging_rgb},
+    io::{bed::read_bed, readers::TabixReader},
     positions::GenomeRange,
 };
 
@@ -206,25 +206,25 @@ pub fn region_fraction_modified(
 ///
 /// Returns an error if the Tabix index cannot be read or a pileup line is
 /// malformed.
-pub fn region_fraction_from_tabix(
-    pileup_bgz: &str,
+fn region_fraction_from_tabix<R: std::io::Read + std::io::Seek>(
+    reader: &mut TabixReader<R>,
     region: &GenomeRange,
 ) -> Result<Option<(f64, u64, u64)>> {
     let mut sum_cov = 0u64;
     let mut sum_mod = 0u64;
 
-    fetch_tabix_lines_with(pileup_bgz, region, |line| {
+    reader.fetch_with(region, |line| {
         let f: Vec<&str> = line.split('\t').collect();
         let nvalid_cov: u64 = f
             .get(9)
-            .ok_or_else(|| anyhow!("Missing Nvalid_cov (col 10) in {pileup_bgz}: {line}"))?
+            .ok_or_else(|| anyhow!("Missing Nvalid_cov (col 10) line : {line}"))?
             .parse()
-            .with_context(|| format!("Bad Nvalid_cov in {pileup_bgz}: {line}"))?;
+            .with_context(|| format!("Bad Nvalid_cov line : {line}"))?;
         let nmod: u64 = f
             .get(11)
-            .ok_or_else(|| anyhow!("Missing Nmod (col 12) in {pileup_bgz}: {line}"))?
+            .ok_or_else(|| anyhow!("Missing Nmod (col 12) line : {line}"))?
             .parse()
-            .with_context(|| format!("Bad Nmod in {pileup_bgz}: {line}"))?;
+            .with_context(|| format!("Bad Nmod in line : {line}"))?;
         sum_cov = sum_cov.saturating_add(nvalid_cov);
         sum_mod = sum_mod.saturating_add(nmod);
         Ok(())
@@ -259,8 +259,8 @@ pub fn region_fraction_from_tabix(
 /// # Errors
 ///
 /// Returns an error if any region has `end <= start` or a Tabix fetch fails.
-pub fn compute_gene_activity(
-    pileup_bgz: &str,
+pub fn compute_gene_activity<R: std::io::Read + std::io::Seek>(
+    reader: &mut TabixReader<R>,
     regions: &[BedRegion],
     min_sum_cov: u64,
 ) -> Result<Vec<GeneActivity>> {
@@ -297,7 +297,7 @@ pub fn compute_gene_activity(
 
         let region = GenomeRange::new(&r.chrom, r.start as u32, r.end as u32);
 
-        let (frac, sum_cov, _) = match region_fraction_from_tabix(pileup_bgz, &region)? {
+        let (frac, sum_cov, _) = match region_fraction_from_tabix(reader, &region)? {
             Some(x) => x,
             None => continue,
         };
@@ -436,8 +436,8 @@ pub fn write_gene_activity_bed9(
 /// # Errors
 ///
 /// Returns an error if any pipeline step fails.
-pub fn pileup_to_activity_bed9(
-    pileup_bgz: &str,
+pub fn pileup_to_activity_bed9<R: std::io::Read + std::io::Seek>(
+    reader: &mut TabixReader<R>,
     regions_bed4_path: &str,
     min_sum_cov: u64,
     out_bed9_path: &str,
@@ -445,7 +445,7 @@ pub fn pileup_to_activity_bed9(
     clip: f64,
 ) -> Result<usize> {
     let regions = read_bed4(regions_bed4_path)?;
-    let activity = compute_gene_activity(pileup_bgz, &regions, min_sum_cov)?;
+    let activity = compute_gene_activity(reader, &regions, min_sum_cov)?;
     let n = activity.len();
     write_gene_activity_bed9(out_bed9_path, track_name, &activity, clip)?;
     Ok(n)
@@ -477,8 +477,9 @@ mod tests {
                 id,
                 config.tumoral_name
             );
+            let mut reader = TabixReader::open(&path)?;
             let n = pileup_to_activity_bed9(
-                &path,
+                &mut reader,
                 regions_bed4_path,
                 100,
                 &out,

src/io/readers.rs

@@ -14,7 +14,7 @@
 use std::{
     fs::{self, File},
     io::{BufRead, BufReader, Read, Write},
-    path::Path,
+    path::Path, sync::Arc,
 };
 
 use anyhow::Context;
@@ -23,6 +23,7 @@ use noodles_bgzf as bgzf;
 use noodles_core::{region::Interval, Position};
 use noodles_csi::BinningIndex;
 use noodles_tabix as tabix;
+use pandora_lib_smb::{SmbReader, SmbSession};
 
 use crate::{
     helpers::TempFileGuard,
@@ -173,111 +174,6 @@ pub fn compress_to_bgzip(input_path: &str) -> anyhow::Result<String> {
     Ok(output_path)
 }
 
-/// Visit every data line in a Tabix-indexed BGZF file that overlaps `region`,
-/// calling `on_line` for each one.
-///
-/// This is the zero-allocation core of the tabix fetch family. A single `String`
-/// buffer is reused across all lines; `on_line` receives a `&str` reference into
-/// that buffer — no intermediate `Vec<String>` is ever built. Callers that need
-/// a collected result (e.g. [`fetch_tabix_lines`], [`fetch_bed`](super::bed::fetch_bed))
-/// simply push into their own accumulator inside the closure.
-///
-/// Header lines (`#`) and blank lines are skipped before `on_line` is called.
-/// Lines are newline-stripped before being passed to the callback.
-///
-/// # Coordinate system
-///
-/// `region` is **0-based half-open** `[start, end)`. The conversion to 1-based
-/// inclusive Tabix positions is performed internally:
-/// - `tabix_start = region.range.start + 1`
-/// - `tabix_end   = region.range.end`
-///
-/// # Arguments
-///
-/// * `bgz_path` - Path to the BGZF file; index expected at `<bgz_path>.tbi`
-/// * `region` - Genomic interval to query (0-based half-open)
-/// * `on_line` - Callback invoked for each matching data line; returning `Err`
-///   aborts the fetch immediately and propagates the error
-///
-/// # Errors
-///
-/// Returns an error if the index or file cannot be read, a seek fails, or
-/// `on_line` returns an error.
-pub fn fetch_tabix_lines_with<F>(
-    bgz_path: &str,
-    region: &GenomeRange,
-    mut on_line: F,
-) -> anyhow::Result<()>
-where
-    F: FnMut(&str) -> anyhow::Result<()>,
-{
-    if region.range.is_empty() {
-        return Ok(());
-    }
-
-    let tbi_path = format!("{bgz_path}.tbi");
-
-    let index = tabix::fs::read(&tbi_path)
-        .with_context(|| format!("cannot read tabix index: {tbi_path}"))?;
-
-    let rname = region.contig();
-
-    let ref_seq_id = match index
-        .header()
-        .ok_or_else(|| anyhow::anyhow!("tabix index has no header: {tbi_path}"))?
-        .reference_sequence_names()
-        .get_index_of(rname.as_bytes())
-    {
-        Some(id) => id,
-        None => return Ok(()),
-    };
-
-    // 0-based half-open [start, end) → 1-based inclusive [start+1, end]
-    let interval_start =
-        Position::try_from(region.range.start as usize + 1).context("region start overflow")?;
-    let interval_end =
-        Position::try_from(region.range.end as usize).context("region end overflow")?;
-
-    let interval = Interval::from(interval_start..=interval_end);
-    let chunks = index
-        .query(ref_seq_id, interval)
-        .context("tabix region query failed")?;
-
-    if chunks.is_empty() {
-        return Ok(());
-    }
-
-    let mut reader = bgzf::io::Reader::new(
-        File::open(bgz_path).with_context(|| format!("cannot open {bgz_path}"))?,
-    );
-
-    // Single buffer reused across all lines — no per-line allocation.
-    let mut buf = String::new();
-
-    for chunk in chunks {
-        reader.seek(chunk.start()).context("BGZF seek failed")?;
-
-        loop {
-            if reader.virtual_position() >= chunk.end() {
-                break;
-            }
-            buf.clear();
-            let n = reader
-                .read_line(&mut buf)
-                .context("BGZF read_line failed")?;
-            if n == 0 {
-                break;
-            }
-            let line = buf.trim_end_matches(|c| c == '\n' || c == '\r');
-            if !line.is_empty() && !line.starts_with('#') {
-                on_line(line)?;
-            }
-        }
-    }
-
-    Ok(())
-}
-
 /// A reusable handle to a Tabix-indexed BGZF file.
 ///
 /// Opening a tabix index and BGZF file is expensive. Use this struct when you
@@ -289,26 +185,68 @@ where
 /// let mut reader = TabixReader::open("/path/to/file.vcf.gz")?;
 /// reader.fetch_with(&region, |line| { println!("{line}"); Ok(()) })?;
 /// ```
-pub struct TabixReader {
+pub struct TabixReader<R> {
     index: tabix::Index,
-    reader: bgzf::io::Reader<File>,
+    reader: bgzf::io::Reader<R>,
     buf: String,
 }
 
-impl TabixReader {
-    /// Open `bgz_path` and its sidecar `.tbi` index.
+impl TabixReader<File> {
+    /// Open `bgz_path` and its sidecar `.tbi` index from local disk.
     pub fn open(bgz_path: &str) -> anyhow::Result<Self> {
         let tbi_path = format!("{bgz_path}.tbi");
         let index = tabix::fs::read(&tbi_path)
             .with_context(|| format!("cannot read tabix index: {tbi_path}"))?;
-        let reader = bgzf::io::Reader::new(
-            File::open(bgz_path).with_context(|| format!("cannot open {bgz_path}"))?,
+        let reader = File::open(bgz_path).with_context(|| format!("cannot open {bgz_path}"))?;
+        Ok(Self::from_reader(reader, index))
+    }
+}
+
+impl TabixReader<SmbReader> {
+    /// Open a tabix-indexed BGZF file hosted over SMB, reusing the session's
+    /// cached SMB reader for the data file and a locally-cached copy of the
+    /// `.tbi` index (parsed via the same on-disk path noodles expects).
+    pub async fn open_smb(session: &Arc<SmbSession>, bgz_rel_path: &str) -> anyhow::Result<Self> {
+        let tbi_rel_path = format!("{bgz_rel_path}.tbi");
+
+        // Download the tabix index bytes and stash them locally — tabix::fs::read
+        // wants a real path, so we materialize one rather than reading in-memory.
+        let tbi_bytes = session
+            .download_file_bytes(&tbi_rel_path)
+            .await
+            .context("download tabix index over SMB")?;
+
+        let cache_dir = session
+            .index_cache_dir()
+            .context("no index cache dir configured for tabix index")?;
+
+        let local_tbi_path = cache_dir.join(
+            bgz_rel_path.replace(['/', '\\'], "_") + ".tbi",
         );
-        Ok(Self {
+
+        fs::write(&local_tbi_path, &tbi_bytes)
+            .with_context(|| format!("write cached tabix index: {}", local_tbi_path.display()))?;
+
+        let index = tabix::fs::read(&local_tbi_path)
+            .with_context(|| format!("parse tabix index: {}", local_tbi_path.display()))?;
+
+        let data = session
+            .open_smb_reader(bgz_rel_path)
+            .await
+            .context("open SMB reader for bgzf data")?;
+
+        Ok(Self::from_reader(data, index))
+    }
+}
+
+impl<R: std::io::Read + std::io::Seek> TabixReader<R> {
+    /// Build a TabixReader from an already-open data reader plus its tabix index bytes.
+    pub fn from_reader(data: R, index: tabix::Index) -> Self {
+        Self {
             index,
-            reader,
+            reader: bgzf::io::Reader::new(data),
             buf: String::new(),
-        })
+        }
     }
 
     /// Visit every data line overlapping `region`, calling `on_line` for each one.
@@ -367,7 +305,7 @@ impl TabixReader {
                 if n == 0 {
                     break;
                 }
-                let line = self.buf.trim_end_matches(|c| c == '\n' || c == '\r');
+                let line = self.buf.trim_end_matches(['\n', '\r']);
                 if !line.is_empty() && !line.starts_with('#') {
                     on_line(line)?;
                 }
@@ -387,9 +325,10 @@ impl TabixReader {
 /// # Errors
 ///
 /// Returns an error if the index or file cannot be read.
-pub fn fetch_tabix_lines(bgz_path: &str, region: &GenomeRange) -> anyhow::Result<Vec<String>> {
+pub fn fetch_tabix_lines<R: std::io::Read + std::io::Seek>(
+    reader: &mut TabixReader<R>, region: &GenomeRange) -> anyhow::Result<Vec<String>> {
     let mut lines = Vec::new();
-    fetch_tabix_lines_with(bgz_path, region, |line| {
+    reader.fetch_with(region, |line| {
         lines.push(line.to_owned());
         Ok(())
     })?;

src/io/vcf.rs

@@ -32,7 +32,7 @@ use crate::{
 
 use super::{
     dict::read_dict,
-    readers::{fetch_tabix_lines_with, get_reader, TabixReader},
+    readers::{get_reader, TabixReader},
 };
 
 /// Load a BGZF-compressed or plain VCF file into memory.
@@ -155,29 +155,29 @@ pub fn write_vcf<B: Borrow<VcfVariant>>(variants: &[B], path: &str) -> anyhow::R
 ///
 /// Returns an error if the index or file cannot be read, or if a fetched line
 /// fails to parse as a VCF variant.
-pub fn fetch_vcf(bgz_path: &str, region: &GenomeRange) -> anyhow::Result<Vec<VcfVariant>> {
-    let mut variants = Vec::new();
-    fetch_tabix_lines_with(bgz_path, region, |line| {
-        let v: VcfVariant = line
-            .parse()
-            .with_context(|| format!("failed to parse VCF line: {line}"))?;
-        if v.position.contig == region.contig && region.range.contains(&v.position.position) {
-            for split in from_multiallelic(v)? {
-                variants.push(split);
-            }
-        }
-        Ok(())
-    })?;
-    Ok(variants)
-}
+// pub fn fetch_vcf(bgz_path: &str, region: &GenomeRange) -> anyhow::Result<Vec<VcfVariant>> {
+//     let mut variants = Vec::new();
+//     fetch_tabix_lines_with(bgz_path, region, |line| {
+//         let v: VcfVariant = line
+//             .parse()
+//             .with_context(|| format!("failed to parse VCF line: {line}"))?;
+//         if v.position.contig == region.contig && region.range.contains(&v.position.position) {
+//             for split in from_multiallelic(v)? {
+//                 variants.push(split);
+//             }
+//         }
+//         Ok(())
+//     })?;
+//     Ok(variants)
+// }
 
 /// Like [`fetch_vcf`] but reuses a pre-opened [`TabixReader`] instead of
 /// reopening the index and BGZF file on every call.
 ///
 /// Use this inside `par_chunks` / `for_each_init` workers where one reader
 /// is initialised per thread and shared across all variants in that chunk.
-pub fn fetch_vcf_with_reader(
-    reader: &mut TabixReader,
+pub fn fetch_vcf_with_reader<R: std::io::Read + std::io::Seek>(
+    reader: &mut TabixReader<R>,
     region: &GenomeRange,
 ) -> anyhow::Result<Vec<VcfVariant>> {
     let mut variants = Vec::new();
@@ -296,7 +296,9 @@ mod tests {
     #[test]
     fn vcf_dbsnp() -> anyhow::Result<()> {
         let config = Config::default();
-        let v = fetch_vcf(&config.db_snp, &GenomeRange::new("chr1", 5656, 5660))?;
+
+        let mut reader = TabixReader::open(&config.db_snp)?;
+        let v = fetch_vcf_with_reader(&mut reader, &GenomeRange::new("chr1", 5656, 5660))?;
         println!("{v:#?}");
         Ok(())
     }

src/lib.rs

@@ -163,7 +163,7 @@ lazy_static! {
 
 #[cfg(test)]
 mod tests {
-    use std::{collections::HashMap, path::Path};
+    use std::{collections::HashMap, path::{Path, PathBuf}};
 
     use annotation::{
         vep::{VepLine, VEP},
@@ -189,7 +189,7 @@ mod tests {
             vcf::VcfCollection,
         },
         helpers::find_files,
-        io::{bed::bedrow_overlaps_par, dict::read_dict, gff::features_ranges},
+        io::{bed::bedrow_overlaps_par, dict::read_dict, fasta::FaiReader, gff::features_ranges},
         positions::{merge_overlapping_genome_ranges, range_intersection_par, sort_ranges},
         scan::scan::somatic_scan,
         variant::{
@@ -996,8 +996,7 @@ mod tests {
         let chr = "chr1";
         let position = 16761;
 
-        let mut fasta_reader =
-            noodles_fasta::io::indexed_reader::Builder::default().build_from_path(c.reference)?;
+        let mut fasta_reader = FaiReader::open(&PathBuf::from(c.reference))?;
         let r = io::fasta::sequence_at(&mut fasta_reader, chr, position, 3)?;
         println!(
             "{r} ({} {:.2})",

src/variant/variant_collection.rs

@@ -9,7 +9,7 @@ use std::{
 use crate::{
     callers::nanomonsv::nanomonsv_insert_classify,
     commands::{exec_jobs, AnyJob},
-    io::{somaticpipe_container::PandoraReader, tsv::TsvLine, vcf::read_vcf},
+    io::{fasta::FaiReader, somaticpipe_container::PandoraReader, tsv::TsvLine, vcf::read_vcf},
     runners::Run,
     variant::vcf_variant::VariantId,
 };
@@ -44,7 +44,7 @@ use crate::{
         TempFileGuard,
     },
     io::{
-        fasta::{open_indexed_fasta, sequence_at},
+        fasta::sequence_at,
         liftover::build_machine_from_chain,
         readers::{get_gz_reader, get_reader},
         vcf::vcf_header,
@@ -382,7 +382,7 @@ impl VariantCollection {
         let chunk_size = self.chunk_size();
 
         self.variants.par_chunks(chunk_size).for_each_init(
-            || noodles_fasta::io::indexed_reader::Builder::default().build_from_path(reference),
+            || FaiReader::open(&PathBuf::from(reference)),
             |reader_res, chunk| {
                 let Ok(ref mut fasta_reader) = reader_res else {
                     error!("Failed to load reference for chunk: {chunk:?}");
@@ -553,8 +553,7 @@ impl VariantCollection {
         self.variants.par_chunks(chunk_size).for_each_init(
             || {
                 let bam = rust_htslib::bam::IndexedReader::from_path(constit_bam_path);
-                let fasta = noodles_fasta::io::indexed_reader::Builder::default()
-                    .build_from_path(&reference);
+                let fasta = FaiReader::open(&PathBuf::from(&reference));
                 (bam, fasta)
             },
             |(bam_res, fasta_res), chunk| {
@@ -722,19 +721,23 @@ impl VariantCollection {
         let n = self.variants.len();
         let mut out: Vec<Option<VcfVariant>> = vec![None; n];
 
-        out.par_chunks_mut(chunk_size).enumerate().try_for_each(
-            |(chunk_idx, out_chunk)| -> anyhow::Result<()> {
-                let mut fasta = open_indexed_fasta(target_fasta_path.as_path())?;
-
-                let start = chunk_idx * chunk_size;
-                let end = (start + out_chunk.len()).min(n);
-
-                for (i, slot) in (start..end).zip(out_chunk.iter_mut()) {
-                    *slot = self.variants[i].liftover_with_fasta(&machine, &mut fasta)?;
-                }
-                Ok(())
-            },
-        )?;
+        // self.variants.par_chunks(chunk_size).for_each_init(
+        //     || FaiReader::open(&PathBuf::from(reference)),
+        //     |reader_res, chunk| {
+        out.par_chunks_mut(chunk_size)
+            .enumerate()
+            .try_for_each_init(
+                || FaiReader::open(&PathBuf::from(target_fasta_path.as_path())),
+                |fa_reader, (chunk_idx, out_chunk)| -> anyhow::Result<()> {
+                    let fasta = fa_reader.as_mut().map_err(|e| anyhow::anyhow!("{e}"))?;
+                    let start = chunk_idx * chunk_size;
+                    let end = (start + out_chunk.len()).min(n);
+                    for (i, slot) in (start..end).zip(out_chunk.iter_mut()) {
+                        *slot = self.variants[i].liftover_with_fasta(&machine, fasta)?;
+                    }
+                    Ok(())
+                },
+            )?;
 
         // Write outputs (serial)
         let mut w_ok = BufWriter::new(File::create(&out_vcf_path)?);
@@ -1654,16 +1657,17 @@ impl Variants {
         let n = self.data.len();
         let mut out: Vec<Option<VcfVariant>> = vec![None; n];
 
-        out.par_chunks_mut(chunk_size).enumerate().try_for_each(
-            |(chunk_idx, out_chunk)| -> anyhow::Result<()> {
-                let mut fasta = open_indexed_fasta(target_fasta_path.as_path())?;
+        out.par_chunks_mut(chunk_size).enumerate().try_for_each_init(
+                || FaiReader::open(&PathBuf::from(target_fasta_path.as_path())),
+                |fa_reader, (chunk_idx, out_chunk)| -> anyhow::Result<()> {
+                let fasta = fa_reader.as_mut().map_err(|e| anyhow::anyhow!("{e}"))?;
 
                 let start = chunk_idx * chunk_size;
                 let end = (start + out_chunk.len()).min(n);
 
                 for (i, slot) in (start..end).zip(out_chunk.iter_mut()) {
                     if let Some(v) = self.data[i].vcf_variants.first() {
-                        *slot = v.liftover_with_fasta(&machine, &mut fasta)?;
+                        *slot = v.liftover_with_fasta(&machine, fasta)?;
                     }
                 }
                 Ok(())

src/variant/vcf_variant.rs

@@ -108,13 +108,12 @@
 
 use crate::{
     annotation::{
-        dbsnp::{DbSnpFreq, DbSnpFreqEntry},
-        Annotations,
+        Annotations, dbsnp::{DbSnpFreq, DbSnpFreqEntry}
     },
-    helpers::{estimate_shannon_entropy, mean, revcomp, Hash128},
-    io::fasta::sequence_range,
+    helpers::{Hash128, estimate_shannon_entropy, mean, revcomp},
+    io::fasta::{FaiReader, sequence_range},
     pipes::ShouldRun,
-    positions::{contig_to_num, GenomePosition, GetGenomePosition, VcfPosition},
+    positions::{GenomePosition, GetGenomePosition, VcfPosition, contig_to_num},
     runners::Run,
     variant::variant_collection::VariantCollection,
 };
@@ -128,7 +127,6 @@ use std::{
     cmp::Ordering,
     collections::{BTreeSet, HashSet},
     fmt,
-    fs::File,
     hash::Hash,
     str::FromStr,
 };
@@ -297,10 +295,10 @@ impl VcfVariant {
     /// - uses REF length to define the lifted interval
     /// - rewrites REF to match target FASTA
     /// - if mapping is on reverse strand, revcomps REF and ALT (sequence alleles)
-    pub fn liftover_with_fasta(
+    pub fn liftover_with_fasta<R: std::io::Read + std::io::Seek>(
         &self,
         machine: &Machine,
-        fasta: &mut noodles_fasta::io::IndexedReader<noodles_fasta::io::BufReader<File>>,
+        fasta: &mut FaiReader<R>,
     ) -> anyhow::Result<Option<Self>> {
         let src_contig = self.position.contig(); // "chr1" etc.
         let src_start0 = self.position.position as u64;