Feature: decompressing for gzip file, supporting flate2 library and shell pipeline approach

Author: Hang Ji

flatgfa/Cargo.toml

@@ -12,6 +12,7 @@ argh = "0.1.13"
 atoi = "2.0.0"
 bit-vec = "0.8.0"
 bstr = "1.12.0"
+flate2 = "1.0"
 memchr = "2.7.4"
 memmap = "0.7.0"
 num_enum = "0.7.3"

flatgfa/src/cli/cmds.rs

@@ -343,13 +343,37 @@ pub fn bed_intersect(args: BEDIntersect) {
 #[derive(FromArgs, PartialEq, Debug)]
 #[argh(subcommand, name = "matrix")]
 pub struct PangenotypeMatrix {
-    /// the GAF file to use for genotyping
-    #[argh(positional)]
-    gaf_file: String,
+    /// the GAF file to use for genotyping (use `-f <file>`)
+    #[argh(option, short = 'f')]
+    file: Option<String>,
+
+    /// read from stdin
+    #[argh(switch)]
+    stdin: bool,
 }
 
 pub fn pangenotype_matrix(gfa: &flatgfa::FlatGFA, args: PangenotypeMatrix) {
-    let matrix = ops::pangenotype::make_pangenotype_matrix(gfa, vec![args.gaf_file]);
+    use std::io::{self, BufRead};
+
+    let matrix = if args.stdin {
+        // Read from stdin stream
+        let stdin = io::stdin();
+        let reader: Box<dyn BufRead> = Box::new(stdin.lock());
+        match ops::pangenotype::make_pangenotype_matrix_from_stream(gfa, reader) {
+            Ok(row) => vec![row],
+            Err(e) => {
+                eprintln!("Error reading from stdin: {}", e);
+                return;
+            }
+        }
+    } else if let Some(file_path) = args.file {
+        // For files (plain or .gz) the existing file-based API handles mmap and gzip
+        ops::pangenotype::make_pangenotype_matrix(gfa, vec![file_path])
+    } else {
+        eprintln!("Please provide a GAF file with -f <file> or use --stdin");
+        return;
+    };
+
     for row in matrix {
         for col in row {
             if col {

flatgfa/src/ops/pangenotype.rs

@@ -1,7 +1,53 @@
 use crate::flatgfa::FlatGFA;
 use crate::memfile;
 use crate::namemap::NameMap;
+use flate2::read::GzDecoder;
 use memchr::memchr;
+use std::fs::File;
+use std::io::{self, BufRead, BufReader, Read};
+
+/// Open a GAF file reader for gzip-compressed files.
+///
+/// This function opens a .gaf.gz file and returns a boxed BufRead trait object
+/// for streaming decompression. Plain text files should use memory mapping instead.
+///
+/// # Arguments
+/// * `path` - The path to a .gaf.gz file
+///
+/// # Returns
+/// A Result containing a boxed trait object implementing BufRead
+fn open_gzip_reader(path: &str) -> io::Result<Box<dyn BufRead>> {
+    let file = File::open(path)?;
+    let decoder = GzDecoder::new(file);
+    Ok(Box::new(BufReader::new(decoder)))
+}
+
+/// Process a GAF stream (implementing BufRead) and update the pangenotype matrix.
+///
+/// This is a generic function that works with any BufRead source (files, gzip streams, stdin, etc.)
+///
+/// # Arguments
+/// * `reader` - A boxed BufRead trait object
+/// * `matrix` - The pangenotype matrix to update
+/// * `file_idx` - The index of this file in the matrix
+/// * `name_map` - The name map for segment lookup
+fn process_gaf_stream(
+    reader: Box<dyn BufRead>,
+    matrix: &mut Vec<Vec<bool>>,
+    file_idx: usize,
+    name_map: &crate::namemap::NameMap,
+) -> io::Result<()> {
+    for line in reader.lines() {
+        let line = line?;
+
+        if line.is_empty() || line.starts_with('#') {
+            continue;
+        }
+
+        process_gaf_line(&line, matrix, file_idx, name_map);
+    }
+    Ok(())
+}
 
 pub fn make_pangenotype_matrix(gfa: &FlatGFA, gaf_files: Vec<String>) -> Vec<Vec<bool>> {
     let num_segments = gfa.segs.len();
@@ -10,56 +56,116 @@ pub fn make_pangenotype_matrix(gfa: &FlatGFA, gaf_files: Vec<String>) -> Vec<Vec
     let name_map = NameMap::build(&gfa);
 
     for (file_idx, gaf_path) in gaf_files.iter().enumerate() {
-        let mmap = memfile::map_file(gaf_path);
-        let mut start = 0;
-
-        while let Some(pos) = memchr(b'\n', &mmap[start..]) {
-            let line_end = start + pos;
-            let line = &mmap[start..line_end];
-            start = line_end + 1;
-
-            if line.is_empty() || line[0] == b'#' {
-                continue;
-            }
-
-            let mut tab_count = 0;
-            let mut idx = 0;
-            while idx < line.len() && tab_count < 5 {
-                if line[idx] == b'\t' {
-                    tab_count += 1;
+        if gaf_path.ends_with(".gz") {
+            // Use BufRead stream for gzip files
+            match open_gzip_reader(gaf_path) {
+                Ok(reader) => {
+                    if let Err(e) = process_gaf_stream(reader, &mut matrix, file_idx, &name_map) {
+                        eprintln!("Error processing GAF stream {}: {}", gaf_path, e);
+                    }
+                }
+                Err(e) => {
+                    eprintln!("Error opening GAF file {}: {}", gaf_path, e);
                 }
-                idx += 1;
             }
+        } else {
+            // Use memory mapping for plain .gaf files (faster)
+            let mmap = memfile::map_file(gaf_path);
+            let mut start = 0;
 
-            if tab_count < 5 || idx >= line.len() {
-                continue;
-            }
-            //The path data is in the 5th field only.
-            let mut end_idx = idx;
-            while end_idx < line.len() && line[end_idx] != b'\t' {
-                end_idx += 1;
-            }
-            let path_field = &line[idx..end_idx];
-
-            // === Parse path field like >12<34>56 ===
-            let mut p = 0;
-            while p < path_field.len() {
-                let byte = path_field[p];
-                if byte == b'>' || byte == b'<' {
-                    p += 1;
-                    let mut num = 0usize;
-                    while p < path_field.len() && path_field[p].is_ascii_digit() {
-                        num = num * 10 + (path_field[p] - b'0') as usize;
-                        p += 1;
-                    }
-                    let seg_id = name_map.get(num);
-                    matrix[file_idx][u32::from(seg_id) as usize] = true;
-                } else {
-                    p += 1;
+            while let Some(pos) = memchr(b'\n', &mmap[start..]) {
+                let line_end = start + pos;
+                let line = &mmap[start..line_end];
+                start = line_end + 1;
+
+                if line.is_empty() || line[0] == b'#' {
+                    continue;
                 }
+
+                process_gaf_line_bytes(line, &mut matrix, file_idx, &name_map);
             }
         }
     }
 
     matrix
 }
+
+pub fn make_pangenotype_matrix_from_stream(
+    gfa: &FlatGFA,
+    reader: Box<dyn BufRead>,
+) -> io::Result<Vec<bool>> {
+    let num_segments = gfa.segs.len();
+    let mut matrix = vec![vec![false; num_segments]];
+    let name_map = NameMap::build(&gfa);
+
+    process_gaf_stream(reader, &mut matrix, 0, &name_map)?;
+
+    Ok(matrix.remove(0))
+}
+
+fn process_gaf_line(line: &str, matrix: &mut Vec<Vec<bool>>, file_idx: usize, name_map: &crate::namemap::NameMap) {
+    let bytes = line.as_bytes();
+    let mut tab_count = 0;
+    let mut idx = 0;
+    while idx < bytes.len() && tab_count < 5 {
+        if bytes[idx] == b'\t' {
+            tab_count += 1;
+        }
+        idx += 1;
+    }
+
+    if tab_count < 5 || idx >= bytes.len() {
+        return;
+    }
+    // The path data is in the 5th field (index 4, after 5 tabs)
+    let mut end_idx = idx;
+    while end_idx < bytes.len() && bytes[end_idx] != b'\t' {
+        end_idx += 1;
+    }
+    let path_field = &bytes[idx..end_idx];
+
+    parse_path_field(path_field, matrix, file_idx, name_map);
+}
+
+fn process_gaf_line_bytes(line: &[u8], matrix: &mut Vec<Vec<bool>>, file_idx: usize, name_map: &crate::namemap::NameMap) {
+    let mut tab_count = 0;
+    let mut idx = 0;
+    while idx < line.len() && tab_count < 5 {
+        if line[idx] == b'\t' {
+            tab_count += 1;
+        }
+        idx += 1;
+    }
+
+    if tab_count < 5 || idx >= line.len() {
+        return;
+    }
+
+    let mut end_idx = idx;
+    while end_idx < line.len() && line[end_idx] != b'\t' {
+        end_idx += 1;
+    }
+    let path_field = &line[idx..end_idx];
+
+    parse_path_field(path_field, matrix, file_idx, name_map);
+}
+
+
+fn parse_path_field(path_field: &[u8], matrix: &mut Vec<Vec<bool>>, file_idx: usize, name_map: &crate::namemap::NameMap) {
+    let mut p = 0;
+    while p < path_field.len() {
+        let byte = path_field[p];
+        if byte == b'>' || byte == b'<' {
+            p += 1;
+            let mut num = 0usize;
+            while p < path_field.len() && path_field[p].is_ascii_digit() {
+                num = num * 10 + (path_field[p] - b'0') as usize;
+                p += 1;
+            }
+            let seg_id = name_map.get(num);
+            matrix[file_idx][u32::from(seg_id) as usize] = true;
+        } else {
+            p += 1;
+        }
+    }
+}