From 724d9798466d6310298dbf7275dfbf493be0fdde Mon Sep 17 00:00:00 2001 From: Hang Ji Date: Wed, 14 May 2025 01:19:02 -0400 Subject: [PATCH 01/14] Adding function to construct pangenotype matrix --- flatgfa-py/flatgfa.pyi | 1 + flatgfa-py/src/lib.rs | 20 ++++++++++++++++++++ 2 files changed, 21 insertions(+) diff --git a/flatgfa-py/flatgfa.pyi b/flatgfa-py/flatgfa.pyi index 16211d43..54b5a009 100644 --- a/flatgfa-py/flatgfa.pyi +++ b/flatgfa-py/flatgfa.pyi @@ -86,6 +86,7 @@ class FlatGFA: def write_gfa(self, filename: str) -> None: ... def all_reads(self, gaf: str) -> GAFParser: ... def print_gaf_lookup(self, gaf: str) -> None: ... + def make_pangenotype_matrix(self,gaf_files: list[str])-> list[list[bool]]:... def parse(filename: str) -> FlatGFA: ... def load(filename: str) -> FlatGFA: ... diff --git a/flatgfa-py/src/lib.rs b/flatgfa-py/src/lib.rs index 446eb6d5..73de71da 100644 --- a/flatgfa-py/src/lib.rs +++ b/flatgfa-py/src/lib.rs @@ -167,6 +167,26 @@ impl PyFlatGFA { println!(); } } + fn make_pangenotype_matrix(&self, gaf_files:Vec) -> Vec>{ + let num_segments = self.segments().__len__(); + //Creation of the matrix, the assumption is the segments are ordered from 0 to 100, etc. + let mut matrix = vec![vec![false; num_segments]; gaf_files.len()]; + let current_gfa = self.0.view(); + let name_map = flatgfa::namemap::NameMap::build(¤t_gfa); + for (file_idx, gaf_file_name) in gaf_files.iter().enumerate(){ + let fake_gaf = self.all_reads(gaf_file_name); + let real_gaf = fake_gaf.gaf_buf.view_gafparser(0); + for lines in real_gaf{ + let path_chunker = flatgfa::ops::gaf::PathChunker::new(¤t_gfa, &name_map, lines); + for event in path_chunker{ + let segment_id = event.handle.segment().index(); + matrix[file_idx][segment_id]= true; + } + } + } + matrix + + } } /// A reference to a list of *any* type within a FlatGFA. From 88966690926c045802d127564057fc03db494594 Mon Sep 17 00:00:00 2001 From: Hang Ji Date: Wed, 14 May 2025 01:34:07 -0400 Subject: [PATCH 02/14] fommating --- flatgfa-py/flatgfa.pyi | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/flatgfa-py/flatgfa.pyi b/flatgfa-py/flatgfa.pyi index 54b5a009..ffbb8eac 100644 --- a/flatgfa-py/flatgfa.pyi +++ b/flatgfa-py/flatgfa.pyi @@ -86,7 +86,7 @@ class FlatGFA: def write_gfa(self, filename: str) -> None: ... def all_reads(self, gaf: str) -> GAFParser: ... def print_gaf_lookup(self, gaf: str) -> None: ... - def make_pangenotype_matrix(self,gaf_files: list[str])-> list[list[bool]]:... + def make_pangenotype_matrix(self, gaf_files: list[str]) -> list[list[bool]]: ... def parse(filename: str) -> FlatGFA: ... def load(filename: str) -> FlatGFA: ... From 0fbd91c5f505af3bc76dc39fc9f1cf88b18a3ade Mon Sep 17 00:00:00 2001 From: Hang Ji Date: Mon, 29 Sep 2025 11:35:37 -0400 Subject: [PATCH 03/14] test --- flatgfa-py/examples/matrix.py | 14 ++++++++++++++ flatgfa-py/examples/matrix_gaf_folder/tiny.gaf | 2 ++ flatgfa-py/examples/matrix_gaf_folder/tiny.gfa | 11 +++++++++++ .../examples/matrix_gaf_folder/tiny_smaller.gaf | 1 + 4 files changed, 28 insertions(+) create mode 100644 flatgfa-py/examples/matrix.py create mode 100644 flatgfa-py/examples/matrix_gaf_folder/tiny.gaf create mode 100644 flatgfa-py/examples/matrix_gaf_folder/tiny.gfa create mode 100644 flatgfa-py/examples/matrix_gaf_folder/tiny_smaller.gaf diff --git a/flatgfa-py/examples/matrix.py b/flatgfa-py/examples/matrix.py new file mode 100644 index 00000000..3ec6d566 --- /dev/null +++ b/flatgfa-py/examples/matrix.py @@ -0,0 +1,14 @@ +import pathlib +import flatgfa +TEST_DIR = pathlib.Path(__file__).parent +# CHR6_GFA = TEST_DIR /"chr6.gfa" +# CHR6_GAF = TEST_DIR /"chr6gaf.gaf" +TEST_GFA = TEST_DIR / "../test/tiny.gfa" +GAF_DIR = (TEST_DIR / "./matrix_gaf_folder").resolve() +graph = flatgfa.parse(str(TEST_GFA)) +gaf = [str(p) for p in GAF_DIR.glob("*.gaf")] +pangenotype_matrix = graph.make_pangenotype_matrix(gaf) +assert len(pangenotype_matrix) == len(gaf) +for gaf_path, row in zip(gaf, pangenotype_matrix): + row01 = [int(b) for b in row] + print(pathlib.Path(gaf_path).name, *row01) \ No newline at end of file diff --git a/flatgfa-py/examples/matrix_gaf_folder/tiny.gaf b/flatgfa-py/examples/matrix_gaf_folder/tiny.gaf new file mode 100644 index 00000000..979fd2ac --- /dev/null +++ b/flatgfa-py/examples/matrix_gaf_folder/tiny.gaf @@ -0,0 +1,2 @@ +foo 12 0 12 + >1>2<4 38 5 17 12 12 0 cg:Z:150M +bar 20 0 20 + >1>2>3 30 7 27 20 20 0 cg:Z:150M diff --git a/flatgfa-py/examples/matrix_gaf_folder/tiny.gfa b/flatgfa-py/examples/matrix_gaf_folder/tiny.gfa new file mode 100644 index 00000000..0a9299a8 --- /dev/null +++ b/flatgfa-py/examples/matrix_gaf_folder/tiny.gfa @@ -0,0 +1,11 @@ +H VN:Z:1.0 +S 1 CAAATAAG +S 2 AAATTTTCTGGAGTTCTAT +S 3 TTG +S 4 CCAACTCTCTG +P one 1+,2+,4- * +P two 1+,2+,3+,4- * +L 1 + 2 + 0M +L 2 + 4 - 0M +L 2 + 3 + 0M +L 3 + 4 - 0M diff --git a/flatgfa-py/examples/matrix_gaf_folder/tiny_smaller.gaf b/flatgfa-py/examples/matrix_gaf_folder/tiny_smaller.gaf new file mode 100644 index 00000000..6b65af5a --- /dev/null +++ b/flatgfa-py/examples/matrix_gaf_folder/tiny_smaller.gaf @@ -0,0 +1 @@ +foo 12 0 12 + >1>2<4 38 5 17 12 12 0 cg:Z:150M \ No newline at end of file From 733c9b9d7bc37d4b77cbcda1fc0533d3726adc0f Mon Sep 17 00:00:00 2001 From: Hang Ji Date: Sun, 16 Nov 2025 20:11:23 -0500 Subject: [PATCH 04/14] Faster building pangenoty matrix and unit test --- flatgfa-py/Cargo.toml | 3 + flatgfa-py/examples/depth.py | 2 +- flatgfa-py/examples/matrix.py | 50 ++++- .../matrix_gaf_folder/tiny_smaller.gaf | 1 - .../examples/{matrix_gaf_folder => }/tiny.gaf | 0 .../examples/{matrix_gaf_folder => }/tiny.gfa | 0 flatgfa-py/src/lib.rs | 173 ++++++++++++++++-- flatgfa-py/test/test_gaf.py | 33 ++++ 8 files changed, 239 insertions(+), 23 deletions(-) delete mode 100644 flatgfa-py/examples/matrix_gaf_folder/tiny_smaller.gaf rename flatgfa-py/examples/{matrix_gaf_folder => }/tiny.gaf (100%) rename flatgfa-py/examples/{matrix_gaf_folder => }/tiny.gfa (100%) diff --git a/flatgfa-py/Cargo.toml b/flatgfa-py/Cargo.toml index 3f5edfd9..ee34503c 100644 --- a/flatgfa-py/Cargo.toml +++ b/flatgfa-py/Cargo.toml @@ -11,3 +11,6 @@ crate-type = ["cdylib"] pyo3 = { version = "0.24.2", features = ["abi3-py38", "multiple-pymethods"] } flatgfa = { path = "../flatgfa" } memmap = "0.7.0" +flate2 = "1.0" +bstr = "1.9" +memchr = "2.7" \ No newline at end of file diff --git a/flatgfa-py/examples/depth.py b/flatgfa-py/examples/depth.py index 74b365ac..27ff6320 100644 --- a/flatgfa-py/examples/depth.py +++ b/flatgfa-py/examples/depth.py @@ -1,7 +1,7 @@ import flatgfa from collections import Counter -graph = flatgfa.parse("../tests/k.gfa") +graph = flatgfa.parse("../test/tiny.gfa") depths = Counter() for path in graph.paths: for step in path: diff --git a/flatgfa-py/examples/matrix.py b/flatgfa-py/examples/matrix.py index 3ec6d566..ad1fa207 100644 --- a/flatgfa-py/examples/matrix.py +++ b/flatgfa-py/examples/matrix.py @@ -1,14 +1,54 @@ import pathlib import flatgfa +from itertools import islice +def iter_bits_ones_dash(path, chunk_bytes: int = 64 * 1024): + ws = b" \t\r\n\v\f" + with open(path, "rb") as f: + while True: + chunk = f.read(chunk_bytes) + if not chunk: + break + for b in chunk: + if b in ws: + continue + if b == ord('1'): + yield True + elif b == ord('0'): + yield False + else: + raise ValueError(f"Invalid character in file: {chr(b)!r} (expected '1' or '0')") + +def compare_row_to_file(row_bits, path) -> tuple[bool, int | None]: + n = len(row_bits) + i = 0 + for fb in iter_bits_ones_dash(path): + if i >= n: + # file has extra bits + raise ValueError("File contains more bits than the row length.") + if fb != bool(row_bits[i]): + return False, i + i += 1 + + if i != n: + # file ended early + raise ValueError(f"File ended early: got {i} bits, expected {n}.") + return True, None + TEST_DIR = pathlib.Path(__file__).parent -# CHR6_GFA = TEST_DIR /"chr6.gfa" -# CHR6_GAF = TEST_DIR /"chr6gaf.gaf" -TEST_GFA = TEST_DIR / "../test/tiny.gfa" +TEST_GFA = TEST_DIR / "./matrix_gaf_folder/Chr3.gfa" GAF_DIR = (TEST_DIR / "./matrix_gaf_folder").resolve() +# MASTER_TXT = TEST_DIR / "./matrix_gaf_folder/Chr3-pangenotype" + graph = flatgfa.parse(str(TEST_GFA)) gaf = [str(p) for p in GAF_DIR.glob("*.gaf")] pangenotype_matrix = graph.make_pangenotype_matrix(gaf) -assert len(pangenotype_matrix) == len(gaf) +# assert len(pangenotype_matrix) == len(gaf) +# first_row = list(pangenotype_matrix[0]) +# ok, where = compare_row_to_file(first_row, str(MASTER_TXT)) +# print("OK (all columns match)" if ok else f"Mismatch at column {where}") +FIRST_N = 100 for gaf_path, row in zip(gaf, pangenotype_matrix): row01 = [int(b) for b in row] - print(pathlib.Path(gaf_path).name, *row01) \ No newline at end of file + # print(pathlib.Path(gaf_path).name, *row01) + first_bits = islice(row, FIRST_N) + print(pathlib.Path(gaf_path).name, *map(int, first_bits)) \ No newline at end of file diff --git a/flatgfa-py/examples/matrix_gaf_folder/tiny_smaller.gaf b/flatgfa-py/examples/matrix_gaf_folder/tiny_smaller.gaf deleted file mode 100644 index 6b65af5a..00000000 --- a/flatgfa-py/examples/matrix_gaf_folder/tiny_smaller.gaf +++ /dev/null @@ -1 +0,0 @@ -foo 12 0 12 + >1>2<4 38 5 17 12 12 0 cg:Z:150M \ No newline at end of file diff --git a/flatgfa-py/examples/matrix_gaf_folder/tiny.gaf b/flatgfa-py/examples/tiny.gaf similarity index 100% rename from flatgfa-py/examples/matrix_gaf_folder/tiny.gaf rename to flatgfa-py/examples/tiny.gaf diff --git a/flatgfa-py/examples/matrix_gaf_folder/tiny.gfa b/flatgfa-py/examples/tiny.gfa similarity index 100% rename from flatgfa-py/examples/matrix_gaf_folder/tiny.gfa rename to flatgfa-py/examples/tiny.gfa diff --git a/flatgfa-py/src/lib.rs b/flatgfa-py/src/lib.rs index 73de71da..b17966bf 100644 --- a/flatgfa-py/src/lib.rs +++ b/flatgfa-py/src/lib.rs @@ -10,6 +10,12 @@ use pyo3::BoundObject; use std::io::Write; use std::str; use std::sync::Arc; +use std::fs::File; +use std::io::{BufRead, BufReader}; +use flate2::read::GzDecoder; +use bstr::{BStr, ByteSlice}; +use memchr::memchr; +use flatgfa::ops::gaf::GAFLine; /// Storage for a FlatGFA. /// @@ -167,26 +173,161 @@ impl PyFlatGFA { println!(); } } - fn make_pangenotype_matrix(&self, gaf_files:Vec) -> Vec>{ - let num_segments = self.segments().__len__(); - //Creation of the matrix, the assumption is the segments are ordered from 0 to 100, etc. - let mut matrix = vec![vec![false; num_segments]; gaf_files.len()]; - let current_gfa = self.0.view(); - let name_map = flatgfa::namemap::NameMap::build(¤t_gfa); - for (file_idx, gaf_file_name) in gaf_files.iter().enumerate(){ - let fake_gaf = self.all_reads(gaf_file_name); - let real_gaf = fake_gaf.gaf_buf.view_gafparser(0); - for lines in real_gaf{ - let path_chunker = flatgfa::ops::gaf::PathChunker::new(¤t_gfa, &name_map, lines); - for event in path_chunker{ - let segment_id = event.handle.segment().index(); - matrix[file_idx][segment_id]= true; + // fn make_pangenotype_matrix(&self, gaf_files:Vec) -> Vec>{ + // let num_segments = self.segments().__len__(); + // //Creation of the matrix, the assumption is the segments are ordered from 0 to 100, etc. + // let mut matrix = vec![vec![false; num_segments]; gaf_files.len()]; + // let current_gfa = self.0.view(); + // let name_map = flatgfa::namemap::NameMap::build(¤t_gfa); + // for (file_idx, gaf_file_name) in gaf_files.iter().enumerate(){ + // let fake_gaf = self.all_reads(gaf_file_name); + // let real_gaf = fake_gaf.gaf_buf.view_gafparser(0); + // for lines in real_gaf{ + // let path_chunker = flatgfa::ops::gaf::PathChunker::new(¤t_gfa, &name_map, lines); + // for event in path_chunker{ + // let segment_id = event.handle.segment().index(); + // matrix[file_idx][segment_id]= true; + // } + // } + // } + // matrix + + // } +// pub fn make_pangenotype_matrix(&self, gaf_files: Vec) -> Vec> { +// let num_segments = self.segments().__len__(); +// let mut matrix = vec![vec![false; num_segments]; gaf_files.len()]; + +// let gfa = self.0.view(); +// let name_map = flatgfa::namemap::NameMap::build(&gfa); + +// for (file_idx, gaf_path) in gaf_files.iter().enumerate() { +// eprintln!("Processing {}", gaf_path); + +// // === 1. Open gzipped file for streaming === +// let file = File::open(gaf_path) +// .unwrap_or_else(|_| panic!("Cannot open GAF file: {}", gaf_path)); +// let decoder = GzDecoder::new(file); + +// // Use a large buffer (1MB) to minimize syscall and decompression overhead +// let mut reader = BufReader::with_capacity(1 << 20, decoder); + +// // === 2. Stream and process in-place === +// loop { +// let buf = reader.fill_buf().unwrap(); +// if buf.is_empty() { +// break; // EOF +// } + +// let mut start = 0; +// while let Some(pos) = memchr(b'\n', &buf[start..]) { +// let line_end = start + pos; +// let line = &buf[start..line_end]; + +// if !line.is_empty() && line[0] != b'#' { +// // === Inline process_line === +// // Split once by tabs, no allocations +// let mut fields = line.split(|&b| b == b'\t'); + +// // Skip first 5 fields +// for _ in 0..5 { +// fields.next(); +// } + +// // 6th field: path +// if let Some(path_field) = fields.next() { +// let gaf_line = flatgfa::ops::gaf::GAFLine { +// name: BStr::new(b"streamed"), +// start: 0, +// end: 0, +// path: path_field, +// }; + +// // Iterate through segments in this path +// let path_chunker = +// flatgfa::ops::gaf::PathChunker::new(&gfa, &name_map, gaf_line); + +// for event in path_chunker { +// let seg_id = event.handle.segment().index(); +// matrix[file_idx][seg_id] = true; +// } +// } +// } + +// start = line_end + 1; // skip past newline +// } + +// // Tell BufReader we processed 'start' bytes +// reader.consume(start); +// } + +// eprintln!("Finished {}", gaf_path); +// } + +// matrix +// } +pub fn make_pangenotype_matrix(&self, gaf_files: Vec) -> Vec> { + let num_segments = self.segments().__len__(); + let mut matrix = vec![vec![false; num_segments]; gaf_files.len()]; + + let gfa = self.0.view(); + let name_map = flatgfa::namemap::NameMap::build(&gfa); + + for (file_idx, gaf_path) in gaf_files.iter().enumerate() { + let file = File::open(gaf_path).expect("Cannot open GAF file"); + let mmap = unsafe { Mmap::map(&file).expect("mmap failed") }; + 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; + } + idx += 1; + } + + 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; } } } - matrix - } + + matrix +} } /// A reference to a list of *any* type within a FlatGFA. diff --git a/flatgfa-py/test/test_gaf.py b/flatgfa-py/test/test_gaf.py index 30d8fb9d..c1c54473 100644 --- a/flatgfa-py/test/test_gaf.py +++ b/flatgfa-py/test/test_gaf.py @@ -1,5 +1,6 @@ import pathlib import flatgfa +import numpy as np TEST_DIR = pathlib.Path(__file__).parent TEST_GFA = TEST_DIR / "tiny.gfa" @@ -24,3 +25,35 @@ def test_gaf_ranges(): [(5, 8), (0, 9), (1, 0)], [(7, 8), (0, 18), (0, 0)], ] + +def test_construct_pangenotype_matrix(): + gfa = flatgfa.parse_bytes(TEST_GFA.read_bytes()) + pangenotype_matrix = gfa.make_pangenotype_matrix([(str(TEST_GAF))]) + assert pangenotype_matrix == [[True, True,True,True]] + +def test_pangenotype_regression(): + gfa = flatgfa.parse_bytes(TEST_GFA.read_bytes()) + pangenotype_matrix = gfa.make_pangenotype_matrix([str(TEST_GAF)]) + + # Optional: inspect structure + print("Matrix:", pangenotype_matrix) + + # Example nested loop: iterate over rows and values + for row in pangenotype_matrix: + for val in row: + print("Segment present:", val) + + # Convert nested list of bools -> NumPy array of floats + a = np.array(pangenotype_matrix, dtype=float) + b = np.random.rand(a.shape[0]) + result = np.linalg.lstsq(a, b, rcond=None) + x, residuals, rank, s = result + + # Print or assert just to verify + print("Coefficients:", x) + print("Residuals:", residuals) + print("Rank:", rank) + print("Singular values:", s) + + # Simple test condition: check the shapes + assert x.shape[0] == a.shape[1] From f1b3649795b04e165dff19b127b9e8b406a331f8 Mon Sep 17 00:00:00 2001 From: Adrian Sampson Date: Sun, 23 Nov 2025 15:10:09 -0500 Subject: [PATCH 05/14] Remove some dead code and reformat We weren't actually using flate2 anywhere yet, so I have removed the dependency. --- Cargo.lock | 2 + flatgfa-py/Cargo.toml | 3 +- flatgfa-py/src/lib.rs | 263 +++++++++++++++++++++--------------------- 3 files changed, 132 insertions(+), 136 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 67d3b6e5..d2add96e 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -129,7 +129,9 @@ dependencies = [ name = "flatgfa-py" version = "0.1.0" dependencies = [ + "bstr", "flatgfa", + "memchr", "memmap", "pyo3", ] diff --git a/flatgfa-py/Cargo.toml b/flatgfa-py/Cargo.toml index ee34503c..32ae365a 100644 --- a/flatgfa-py/Cargo.toml +++ b/flatgfa-py/Cargo.toml @@ -11,6 +11,5 @@ crate-type = ["cdylib"] pyo3 = { version = "0.24.2", features = ["abi3-py38", "multiple-pymethods"] } flatgfa = { path = "../flatgfa" } memmap = "0.7.0" -flate2 = "1.0" bstr = "1.9" -memchr = "2.7" \ No newline at end of file +memchr = "2.7" diff --git a/flatgfa-py/src/lib.rs b/flatgfa-py/src/lib.rs index b17966bf..ef20ff76 100644 --- a/flatgfa-py/src/lib.rs +++ b/flatgfa-py/src/lib.rs @@ -2,20 +2,16 @@ use flatgfa::namemap::NameMap; use flatgfa::ops::gaf::{ChunkEvent, GAFParser}; use flatgfa::pool::Id; use flatgfa::{self, file, memfile, print, FlatGFA, Handle, HeapGFAStore}; +use memchr::memchr; use memmap::Mmap; use pyo3::exceptions::PyIndexError; use pyo3::prelude::*; use pyo3::types::{PyBytes, PySlice}; use pyo3::BoundObject; +use std::fs::File; use std::io::Write; use std::str; use std::sync::Arc; -use std::fs::File; -use std::io::{BufRead, BufReader}; -use flate2::read::GzDecoder; -use bstr::{BStr, ByteSlice}; -use memchr::memchr; -use flatgfa::ops::gaf::GAFLine; /// Storage for a FlatGFA. /// @@ -175,7 +171,7 @@ impl PyFlatGFA { } // fn make_pangenotype_matrix(&self, gaf_files:Vec) -> Vec>{ // let num_segments = self.segments().__len__(); - // //Creation of the matrix, the assumption is the segments are ordered from 0 to 100, etc. + // //Creation of the matrix, the assumption is the segments are ordered from 0 to 100, etc. // let mut matrix = vec![vec![false; num_segments]; gaf_files.len()]; // let current_gfa = self.0.view(); // let name_map = flatgfa::namemap::NameMap::build(¤t_gfa); @@ -191,143 +187,142 @@ impl PyFlatGFA { // } // } // matrix - + // } -// pub fn make_pangenotype_matrix(&self, gaf_files: Vec) -> Vec> { -// let num_segments = self.segments().__len__(); -// let mut matrix = vec![vec![false; num_segments]; gaf_files.len()]; - -// let gfa = self.0.view(); -// let name_map = flatgfa::namemap::NameMap::build(&gfa); - -// for (file_idx, gaf_path) in gaf_files.iter().enumerate() { -// eprintln!("Processing {}", gaf_path); - -// // === 1. Open gzipped file for streaming === -// let file = File::open(gaf_path) -// .unwrap_or_else(|_| panic!("Cannot open GAF file: {}", gaf_path)); -// let decoder = GzDecoder::new(file); - -// // Use a large buffer (1MB) to minimize syscall and decompression overhead -// let mut reader = BufReader::with_capacity(1 << 20, decoder); - -// // === 2. Stream and process in-place === -// loop { -// let buf = reader.fill_buf().unwrap(); -// if buf.is_empty() { -// break; // EOF -// } - -// let mut start = 0; -// while let Some(pos) = memchr(b'\n', &buf[start..]) { -// let line_end = start + pos; -// let line = &buf[start..line_end]; - -// if !line.is_empty() && line[0] != b'#' { -// // === Inline process_line === -// // Split once by tabs, no allocations -// let mut fields = line.split(|&b| b == b'\t'); - -// // Skip first 5 fields -// for _ in 0..5 { -// fields.next(); -// } - -// // 6th field: path -// if let Some(path_field) = fields.next() { -// let gaf_line = flatgfa::ops::gaf::GAFLine { -// name: BStr::new(b"streamed"), -// start: 0, -// end: 0, -// path: path_field, -// }; - -// // Iterate through segments in this path -// let path_chunker = -// flatgfa::ops::gaf::PathChunker::new(&gfa, &name_map, gaf_line); - -// for event in path_chunker { -// let seg_id = event.handle.segment().index(); -// matrix[file_idx][seg_id] = true; -// } -// } -// } - -// start = line_end + 1; // skip past newline -// } - -// // Tell BufReader we processed 'start' bytes -// reader.consume(start); -// } - -// eprintln!("Finished {}", gaf_path); -// } - -// matrix -// } -pub fn make_pangenotype_matrix(&self, gaf_files: Vec) -> Vec> { - let num_segments = self.segments().__len__(); - let mut matrix = vec![vec![false; num_segments]; gaf_files.len()]; - - let gfa = self.0.view(); - let name_map = flatgfa::namemap::NameMap::build(&gfa); - - for (file_idx, gaf_path) in gaf_files.iter().enumerate() { - let file = File::open(gaf_path).expect("Cannot open GAF file"); - let mmap = unsafe { Mmap::map(&file).expect("mmap failed") }; - 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; - } + // pub fn make_pangenotype_matrix(&self, gaf_files: Vec) -> Vec> { + // let num_segments = self.segments().__len__(); + // let mut matrix = vec![vec![false; num_segments]; gaf_files.len()]; + + // let gfa = self.0.view(); + // let name_map = flatgfa::namemap::NameMap::build(&gfa); + + // for (file_idx, gaf_path) in gaf_files.iter().enumerate() { + // eprintln!("Processing {}", gaf_path); + + // // === 1. Open gzipped file for streaming === + // let file = File::open(gaf_path) + // .unwrap_or_else(|_| panic!("Cannot open GAF file: {}", gaf_path)); + // let decoder = GzDecoder::new(file); + + // // Use a large buffer (1MB) to minimize syscall and decompression overhead + // let mut reader = BufReader::with_capacity(1 << 20, decoder); + + // // === 2. Stream and process in-place === + // loop { + // let buf = reader.fill_buf().unwrap(); + // if buf.is_empty() { + // break; // EOF + // } + + // let mut start = 0; + // while let Some(pos) = memchr(b'\n', &buf[start..]) { + // let line_end = start + pos; + // let line = &buf[start..line_end]; + + // if !line.is_empty() && line[0] != b'#' { + // // === Inline process_line === + // // Split once by tabs, no allocations + // let mut fields = line.split(|&b| b == b'\t'); + + // // Skip first 5 fields + // for _ in 0..5 { + // fields.next(); + // } + + // // 6th field: path + // if let Some(path_field) = fields.next() { + // let gaf_line = flatgfa::ops::gaf::GAFLine { + // name: BStr::new(b"streamed"), + // start: 0, + // end: 0, + // path: path_field, + // }; + + // // Iterate through segments in this path + // let path_chunker = + // flatgfa::ops::gaf::PathChunker::new(&gfa, &name_map, gaf_line); + + // for event in path_chunker { + // let seg_id = event.handle.segment().index(); + // matrix[file_idx][seg_id] = true; + // } + // } + // } + + // start = line_end + 1; // skip past newline + // } + + // // Tell BufReader we processed 'start' bytes + // reader.consume(start); + // } + + // eprintln!("Finished {}", gaf_path); + // } + + // matrix + // } + pub fn make_pangenotype_matrix(&self, gaf_files: Vec) -> Vec> { + let num_segments = self.segments().__len__(); + let mut matrix = vec![vec![false; num_segments]; gaf_files.len()]; + + let gfa = self.0.view(); + let name_map = flatgfa::namemap::NameMap::build(&gfa); + + for (file_idx, gaf_path) in gaf_files.iter().enumerate() { + let file = File::open(gaf_path).expect("Cannot open GAF file"); + let mmap = unsafe { Mmap::map(&file).expect("mmap failed") }; + 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; - 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 line.is_empty() || line[0] == b'#' { + continue; } - idx += 1; - } - 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 mut tab_count = 0; + let mut idx = 0; + while idx < line.len() && tab_count < 5 { + if line[idx] == b'\t' { + tab_count += 1; } - let seg_id = name_map.get(num); + idx += 1; + } + + 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; + } else { + p += 1; + } } } } - } - matrix -} + matrix + } } /// A reference to a list of *any* type within a FlatGFA. From 93fe0ba26c406765b9efa44d8dafecea1bc019d1 Mon Sep 17 00:00:00 2001 From: Adrian Sampson Date: Sun, 23 Nov 2025 15:11:12 -0500 Subject: [PATCH 06/14] Remove one more dead dependency --- Cargo.lock | 1 - flatgfa-py/Cargo.toml | 1 - 2 files changed, 2 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index d2add96e..2184758d 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -129,7 +129,6 @@ dependencies = [ name = "flatgfa-py" version = "0.1.0" dependencies = [ - "bstr", "flatgfa", "memchr", "memmap", diff --git a/flatgfa-py/Cargo.toml b/flatgfa-py/Cargo.toml index 32ae365a..7364baa7 100644 --- a/flatgfa-py/Cargo.toml +++ b/flatgfa-py/Cargo.toml @@ -11,5 +11,4 @@ crate-type = ["cdylib"] pyo3 = { version = "0.24.2", features = ["abi3-py38", "multiple-pymethods"] } flatgfa = { path = "../flatgfa" } memmap = "0.7.0" -bstr = "1.9" memchr = "2.7" From 6f598d9f6906273a1bf17cde81431d0ccaef223c Mon Sep 17 00:00:00 2001 From: Adrian Sampson Date: Sun, 23 Nov 2025 15:11:41 -0500 Subject: [PATCH 07/14] Remove commented-out dead code Trying to keep the focus on the current version... --- flatgfa-py/src/lib.rs | 92 ------------------------------------------- 1 file changed, 92 deletions(-) diff --git a/flatgfa-py/src/lib.rs b/flatgfa-py/src/lib.rs index ef20ff76..513e3ca0 100644 --- a/flatgfa-py/src/lib.rs +++ b/flatgfa-py/src/lib.rs @@ -169,98 +169,6 @@ impl PyFlatGFA { println!(); } } - // fn make_pangenotype_matrix(&self, gaf_files:Vec) -> Vec>{ - // let num_segments = self.segments().__len__(); - // //Creation of the matrix, the assumption is the segments are ordered from 0 to 100, etc. - // let mut matrix = vec![vec![false; num_segments]; gaf_files.len()]; - // let current_gfa = self.0.view(); - // let name_map = flatgfa::namemap::NameMap::build(¤t_gfa); - // for (file_idx, gaf_file_name) in gaf_files.iter().enumerate(){ - // let fake_gaf = self.all_reads(gaf_file_name); - // let real_gaf = fake_gaf.gaf_buf.view_gafparser(0); - // for lines in real_gaf{ - // let path_chunker = flatgfa::ops::gaf::PathChunker::new(¤t_gfa, &name_map, lines); - // for event in path_chunker{ - // let segment_id = event.handle.segment().index(); - // matrix[file_idx][segment_id]= true; - // } - // } - // } - // matrix - - // } - // pub fn make_pangenotype_matrix(&self, gaf_files: Vec) -> Vec> { - // let num_segments = self.segments().__len__(); - // let mut matrix = vec![vec![false; num_segments]; gaf_files.len()]; - - // let gfa = self.0.view(); - // let name_map = flatgfa::namemap::NameMap::build(&gfa); - - // for (file_idx, gaf_path) in gaf_files.iter().enumerate() { - // eprintln!("Processing {}", gaf_path); - - // // === 1. Open gzipped file for streaming === - // let file = File::open(gaf_path) - // .unwrap_or_else(|_| panic!("Cannot open GAF file: {}", gaf_path)); - // let decoder = GzDecoder::new(file); - - // // Use a large buffer (1MB) to minimize syscall and decompression overhead - // let mut reader = BufReader::with_capacity(1 << 20, decoder); - - // // === 2. Stream and process in-place === - // loop { - // let buf = reader.fill_buf().unwrap(); - // if buf.is_empty() { - // break; // EOF - // } - - // let mut start = 0; - // while let Some(pos) = memchr(b'\n', &buf[start..]) { - // let line_end = start + pos; - // let line = &buf[start..line_end]; - - // if !line.is_empty() && line[0] != b'#' { - // // === Inline process_line === - // // Split once by tabs, no allocations - // let mut fields = line.split(|&b| b == b'\t'); - - // // Skip first 5 fields - // for _ in 0..5 { - // fields.next(); - // } - - // // 6th field: path - // if let Some(path_field) = fields.next() { - // let gaf_line = flatgfa::ops::gaf::GAFLine { - // name: BStr::new(b"streamed"), - // start: 0, - // end: 0, - // path: path_field, - // }; - - // // Iterate through segments in this path - // let path_chunker = - // flatgfa::ops::gaf::PathChunker::new(&gfa, &name_map, gaf_line); - - // for event in path_chunker { - // let seg_id = event.handle.segment().index(); - // matrix[file_idx][seg_id] = true; - // } - // } - // } - - // start = line_end + 1; // skip past newline - // } - - // // Tell BufReader we processed 'start' bytes - // reader.consume(start); - // } - - // eprintln!("Finished {}", gaf_path); - // } - - // matrix - // } pub fn make_pangenotype_matrix(&self, gaf_files: Vec) -> Vec> { let num_segments = self.segments().__len__(); let mut matrix = vec![vec![false; num_segments]; gaf_files.len()]; From baa4123274db956613a45278d2069f114755bf4c Mon Sep 17 00:00:00 2001 From: Adrian Sampson Date: Sun, 23 Nov 2025 15:15:56 -0500 Subject: [PATCH 08/14] Remove an unrelated change --- flatgfa-py/examples/depth.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/flatgfa-py/examples/depth.py b/flatgfa-py/examples/depth.py index 27ff6320..74b365ac 100644 --- a/flatgfa-py/examples/depth.py +++ b/flatgfa-py/examples/depth.py @@ -1,7 +1,7 @@ import flatgfa from collections import Counter -graph = flatgfa.parse("../test/tiny.gfa") +graph = flatgfa.parse("../tests/k.gfa") depths = Counter() for path in graph.paths: for step in path: From bf765c8aa46bb50fdacdde863080f9aec534205c Mon Sep 17 00:00:00 2001 From: Adrian Sampson Date: Sun, 23 Nov 2025 15:16:24 -0500 Subject: [PATCH 09/14] Reformat Python files --- flatgfa-py/examples/matrix.py | 19 +++++++++++++------ flatgfa-py/test/test_gaf.py | 8 ++++++-- 2 files changed, 19 insertions(+), 8 deletions(-) diff --git a/flatgfa-py/examples/matrix.py b/flatgfa-py/examples/matrix.py index ad1fa207..91ae50e6 100644 --- a/flatgfa-py/examples/matrix.py +++ b/flatgfa-py/examples/matrix.py @@ -1,6 +1,9 @@ import pathlib -import flatgfa from itertools import islice + +import flatgfa + + def iter_bits_ones_dash(path, chunk_bytes: int = 64 * 1024): ws = b" \t\r\n\v\f" with open(path, "rb") as f: @@ -11,12 +14,15 @@ def iter_bits_ones_dash(path, chunk_bytes: int = 64 * 1024): for b in chunk: if b in ws: continue - if b == ord('1'): + if b == ord("1"): yield True - elif b == ord('0'): + elif b == ord("0"): yield False else: - raise ValueError(f"Invalid character in file: {chr(b)!r} (expected '1' or '0')") + raise ValueError( + f"Invalid character in file: {chr(b)!r} (expected '1' or '0')" + ) + def compare_row_to_file(row_bits, path) -> tuple[bool, int | None]: n = len(row_bits) @@ -34,6 +40,7 @@ def compare_row_to_file(row_bits, path) -> tuple[bool, int | None]: raise ValueError(f"File ended early: got {i} bits, expected {n}.") return True, None + TEST_DIR = pathlib.Path(__file__).parent TEST_GFA = TEST_DIR / "./matrix_gaf_folder/Chr3.gfa" GAF_DIR = (TEST_DIR / "./matrix_gaf_folder").resolve() @@ -48,7 +55,7 @@ def compare_row_to_file(row_bits, path) -> tuple[bool, int | None]: # print("OK (all columns match)" if ok else f"Mismatch at column {where}") FIRST_N = 100 for gaf_path, row in zip(gaf, pangenotype_matrix): - row01 = [int(b) for b in row] + row01 = [int(b) for b in row] # print(pathlib.Path(gaf_path).name, *row01) first_bits = islice(row, FIRST_N) - print(pathlib.Path(gaf_path).name, *map(int, first_bits)) \ No newline at end of file + print(pathlib.Path(gaf_path).name, *map(int, first_bits)) diff --git a/flatgfa-py/test/test_gaf.py b/flatgfa-py/test/test_gaf.py index c1c54473..55a135bd 100644 --- a/flatgfa-py/test/test_gaf.py +++ b/flatgfa-py/test/test_gaf.py @@ -1,7 +1,9 @@ import pathlib -import flatgfa + import numpy as np +import flatgfa + TEST_DIR = pathlib.Path(__file__).parent TEST_GFA = TEST_DIR / "tiny.gfa" TEST_GAF = TEST_DIR / "tiny.gaf" @@ -26,10 +28,12 @@ def test_gaf_ranges(): [(7, 8), (0, 18), (0, 0)], ] + def test_construct_pangenotype_matrix(): gfa = flatgfa.parse_bytes(TEST_GFA.read_bytes()) pangenotype_matrix = gfa.make_pangenotype_matrix([(str(TEST_GAF))]) - assert pangenotype_matrix == [[True, True,True,True]] + assert pangenotype_matrix == [[True, True, True, True]] + def test_pangenotype_regression(): gfa = flatgfa.parse_bytes(TEST_GFA.read_bytes()) From a51c5e28cd912783faab720297dd9afcd25ae637 Mon Sep 17 00:00:00 2001 From: Adrian Sampson Date: Sun, 23 Nov 2025 15:26:29 -0500 Subject: [PATCH 10/14] Move pangenotype matrix construction to library Instead of putting this directly in the Python bindings, now this is available in the core library. This way, we can also add a command-line way to do the same thing, for instance. --- Cargo.lock | 1 - flatgfa-py/Cargo.toml | 1 - flatgfa-py/src/lib.rs | 64 ++------------------------------- flatgfa/src/ops/mod.rs | 1 + flatgfa/src/ops/pangenotype.rs | 65 ++++++++++++++++++++++++++++++++++ 5 files changed, 69 insertions(+), 63 deletions(-) create mode 100644 flatgfa/src/ops/pangenotype.rs diff --git a/Cargo.lock b/Cargo.lock index 2184758d..67d3b6e5 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -130,7 +130,6 @@ name = "flatgfa-py" version = "0.1.0" dependencies = [ "flatgfa", - "memchr", "memmap", "pyo3", ] diff --git a/flatgfa-py/Cargo.toml b/flatgfa-py/Cargo.toml index 7364baa7..3f5edfd9 100644 --- a/flatgfa-py/Cargo.toml +++ b/flatgfa-py/Cargo.toml @@ -11,4 +11,3 @@ crate-type = ["cdylib"] pyo3 = { version = "0.24.2", features = ["abi3-py38", "multiple-pymethods"] } flatgfa = { path = "../flatgfa" } memmap = "0.7.0" -memchr = "2.7" diff --git a/flatgfa-py/src/lib.rs b/flatgfa-py/src/lib.rs index 513e3ca0..68a85414 100644 --- a/flatgfa-py/src/lib.rs +++ b/flatgfa-py/src/lib.rs @@ -1,14 +1,13 @@ use flatgfa::namemap::NameMap; use flatgfa::ops::gaf::{ChunkEvent, GAFParser}; +use flatgfa::ops::pangenotype; use flatgfa::pool::Id; use flatgfa::{self, file, memfile, print, FlatGFA, Handle, HeapGFAStore}; -use memchr::memchr; use memmap::Mmap; use pyo3::exceptions::PyIndexError; use pyo3::prelude::*; use pyo3::types::{PyBytes, PySlice}; use pyo3::BoundObject; -use std::fs::File; use std::io::Write; use std::str; use std::sync::Arc; @@ -169,67 +168,10 @@ impl PyFlatGFA { println!(); } } - pub fn make_pangenotype_matrix(&self, gaf_files: Vec) -> Vec> { - let num_segments = self.segments().__len__(); - let mut matrix = vec![vec![false; num_segments]; gaf_files.len()]; + pub fn make_pangenotype_matrix(&self, gaf_files: Vec) -> Vec> { let gfa = self.0.view(); - let name_map = flatgfa::namemap::NameMap::build(&gfa); - - for (file_idx, gaf_path) in gaf_files.iter().enumerate() { - let file = File::open(gaf_path).expect("Cannot open GAF file"); - let mmap = unsafe { Mmap::map(&file).expect("mmap failed") }; - 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; - } - idx += 1; - } - - 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; - } - } - } - } - - matrix + pangenotype::make_pangenotype_matrix(&gfa, gaf_files) } } diff --git a/flatgfa/src/ops/mod.rs b/flatgfa/src/ops/mod.rs index 5eaaf835..e6d4773f 100644 --- a/flatgfa/src/ops/mod.rs +++ b/flatgfa/src/ops/mod.rs @@ -3,4 +3,5 @@ pub mod chop; pub mod depth; pub mod extract; pub mod gaf; +pub mod pangenotype; pub mod position; diff --git a/flatgfa/src/ops/pangenotype.rs b/flatgfa/src/ops/pangenotype.rs new file mode 100644 index 00000000..9ff87629 --- /dev/null +++ b/flatgfa/src/ops/pangenotype.rs @@ -0,0 +1,65 @@ +use crate::flatgfa::FlatGFA; +use crate::memfile; +use crate::namemap::NameMap; +use memchr::memchr; + +pub fn make_pangenotype_matrix(gfa: &FlatGFA, gaf_files: Vec) -> Vec> { + let num_segments = gfa.segs.len(); + let mut matrix = vec![vec![false; num_segments]; gaf_files.len()]; + + 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; + } + idx += 1; + } + + 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; + } + } + } + } + + matrix +} From cc48cabfc0bac3222d11b71abbf28af5dbd19723 Mon Sep 17 00:00:00 2001 From: Adrian Sampson Date: Sun, 23 Nov 2025 15:35:29 -0500 Subject: [PATCH 11/14] Add CLI entry point for pangenotype matrix --- flatgfa/src/cli/cmds.rs | 23 +++++++++++++++++++++++ flatgfa/src/cli/main.rs | 4 ++++ 2 files changed, 27 insertions(+) diff --git a/flatgfa/src/cli/cmds.rs b/flatgfa/src/cli/cmds.rs index 61f8bb54..a4bcf449 100644 --- a/flatgfa/src/cli/cmds.rs +++ b/flatgfa/src/cli/cmds.rs @@ -338,3 +338,26 @@ pub fn bed_intersect(args: BEDIntersect) { } } } + +/// construct a pangenotype matrix +#[derive(FromArgs, PartialEq, Debug)] +#[argh(subcommand, name = "matrix")] +pub struct PangenotypeMatrix { + /// the GAF file to use for genotyping + #[argh(positional)] + gaf_file: String, +} + +pub fn pangenotype_matrix(gfa: &flatgfa::FlatGFA, args: PangenotypeMatrix) { + let matrix = ops::pangenotype::make_pangenotype_matrix(gfa, vec![args.gaf_file]); + for row in matrix { + for col in row { + if col { + print!("1"); + } else { + print!("0"); + } + } + println!(); + } +} diff --git a/flatgfa/src/cli/main.rs b/flatgfa/src/cli/main.rs index c2bb0715..c39bebd5 100644 --- a/flatgfa/src/cli/main.rs +++ b/flatgfa/src/cli/main.rs @@ -44,6 +44,7 @@ enum Command { GafLookup(cmds::GAFLookup), Bench(cmds::Bench), BedIntersect(cmds::BEDIntersect), + PangenotypeMatrix(cmds::PangenotypeMatrix), } fn main() -> Result<(), &'static str> { @@ -146,6 +147,9 @@ fn main() -> Result<(), &'static str> { Some(Command::BedIntersect(_sub_args)) => { panic!("Unreachable code"); } + Some(Command::PangenotypeMatrix(sub_args)) => { + cmds::pangenotype_matrix(&gfa, sub_args); + } None => { // Just emit the GFA or FlatGFA file. dump(&gfa, &args.output); From 7dd33141474a740704053588ca377896a6763797 Mon Sep 17 00:00:00 2001 From: Hang Ji Date: Mon, 1 Dec 2025 00:10:33 -0500 Subject: [PATCH 12/14] Feature: decompressing for gzip file, supporting flate2 library and shell pipeline approach --- flatgfa/Cargo.toml | 1 + flatgfa/src/cli/cmds.rs | 32 +++++- flatgfa/src/ops/pangenotype.rs | 192 +++++++++++++++++++++++++-------- 3 files changed, 178 insertions(+), 47 deletions(-) diff --git a/flatgfa/Cargo.toml b/flatgfa/Cargo.toml index faa8d978..3141259a 100644 --- a/flatgfa/Cargo.toml +++ b/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" diff --git a/flatgfa/src/cli/cmds.rs b/flatgfa/src/cli/cmds.rs index a4bcf449..39952adc 100644 --- a/flatgfa/src/cli/cmds.rs +++ b/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 `) + #[argh(option, short = 'f')] + file: Option, + + /// 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 = 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 or use --stdin"); + return; + }; + for row in matrix { for col in row { if col { diff --git a/flatgfa/src/ops/pangenotype.rs b/flatgfa/src/ops/pangenotype.rs index 9ff87629..057e34b9 100644 --- a/flatgfa/src/ops/pangenotype.rs +++ b/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> { + 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, + matrix: &mut Vec>, + 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) -> Vec> { let num_segments = gfa.segs.len(); @@ -10,56 +56,116 @@ pub fn make_pangenotype_matrix(gfa: &FlatGFA, gaf_files: Vec) -> Vec { + 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, +) -> io::Result> { + 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>, 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>, 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>, 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; + } + } +} From 3f503b7ac12f32252e879cf6e80cb51d498ac39d Mon Sep 17 00:00:00 2001 From: Hang Ji Date: Sun, 14 Dec 2025 01:25:26 -0500 Subject: [PATCH 13/14] formatting --- flatgfa/src/ops/pangenotype.rs | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/flatgfa/src/ops/pangenotype.rs b/flatgfa/src/ops/pangenotype.rs index 057e34b9..5cdec7e6 100644 --- a/flatgfa/src/ops/pangenotype.rs +++ b/flatgfa/src/ops/pangenotype.rs @@ -103,7 +103,12 @@ pub fn make_pangenotype_matrix_from_stream( Ok(matrix.remove(0)) } -fn process_gaf_line(line: &str, matrix: &mut Vec>, file_idx: usize, name_map: &crate::namemap::NameMap) { +fn process_gaf_line( + line: &str, + matrix: &mut Vec>, + file_idx: usize, + name_map: &crate::namemap::NameMap, +) { let bytes = line.as_bytes(); let mut tab_count = 0; let mut idx = 0; @@ -127,7 +132,12 @@ fn process_gaf_line(line: &str, matrix: &mut Vec>, file_idx: usize, na parse_path_field(path_field, matrix, file_idx, name_map); } -fn process_gaf_line_bytes(line: &[u8], matrix: &mut Vec>, file_idx: usize, name_map: &crate::namemap::NameMap) { +fn process_gaf_line_bytes( + line: &[u8], + matrix: &mut Vec>, + file_idx: usize, + name_map: &crate::namemap::NameMap, +) { let mut tab_count = 0; let mut idx = 0; while idx < line.len() && tab_count < 5 { @@ -150,8 +160,12 @@ fn process_gaf_line_bytes(line: &[u8], matrix: &mut Vec>, file_idx: us parse_path_field(path_field, matrix, file_idx, name_map); } - -fn parse_path_field(path_field: &[u8], matrix: &mut Vec>, file_idx: usize, name_map: &crate::namemap::NameMap) { +fn parse_path_field( + path_field: &[u8], + matrix: &mut Vec>, + file_idx: usize, + name_map: &crate::namemap::NameMap, +) { let mut p = 0; while p < path_field.len() { let byte = path_field[p]; From bb91e67ccf8174e2569b261a3a70b52d74faafc5 Mon Sep 17 00:00:00 2001 From: Hang Ji Date: Sun, 14 Dec 2025 01:36:43 -0500 Subject: [PATCH 14/14] format --- flatgfa/src/ops/pangenotype.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/flatgfa/src/ops/pangenotype.rs b/flatgfa/src/ops/pangenotype.rs index 5cdec7e6..97c9893d 100644 --- a/flatgfa/src/ops/pangenotype.rs +++ b/flatgfa/src/ops/pangenotype.rs @@ -4,7 +4,7 @@ use crate::namemap::NameMap; use flate2::read::GzDecoder; use memchr::memchr; use std::fs::File; -use std::io::{self, BufRead, BufReader, Read}; +use std::io::{self, BufRead, BufReader}; /// Open a GAF file reader for gzip-compressed files. ///