Merge pull request #99 from pangenome/impg_lace

Add `impg lace` command

Author: AndreaGuarracino

Cargo.lock

@@ -12,7 +12,7 @@ clap = { version = "4.5.42", features = ["derive"] }
 coitrees = "0.4.0"
 num_cpus = "1.17.0"
 rayon = "1.10.0"
-noodles = { version = "0.99.0", features = ["bgzf"] }
+noodles = { version = "0.100.0", features = ["bgzf"] }
 regex = "1.11.1"
 log = "0.4.27"
 env_logger = "0.11.8"
@@ -28,6 +28,14 @@ nalgebra = "0.33.2"
 
 agc-rs = { git = "https://github.com/pangenome/agc-rs.git", optional = true }
 
+# Additional dependencies for lace functionality
+handlegraph = { git = "https://github.com/chfi/rs-handlegraph", rev = "3ac575e4216ce16a16667503a8875e469a40a97a" }
+bitvec = "1.0.1"
+niffler = "2.7.0"
+tempfile = "3.20.0"
+zstd = { version = "0.13.3", features = ["zstdmt"] }
+gzp = "1.0.1"
+
 [features]
 default = ["agc"]
 agc = ["agc-rs"]

Cargo.toml

@@ -183,6 +183,71 @@ Print alignment statistics:
 impg stats -p alignments.paf
 ```
 
+### Lace
+
+Combine multiple GFA files into a unified pangenome graph:
+
+```bash
+# Combine multiple GFA files (you can mix compressed and uncompressed files)
+impg lace -g file1.gfa file2.gfa file3.gfa -o combined.gfa
+
+# Use a list file containing GFA paths
+impg lace -l gfa_files.txt -o combined.gfa
+
+# Fill gaps between contiguous path segments
+impg lace -g *.gfa -o combined.gfa --fill-gaps 1 # Fill with N's
+impg lace -g *.gfa -o combined.gfa --fill-gaps 1 --sequence-files sequence.fa # Fill with sequences
+
+# Fill all gaps, including start and end gaps (end gap filling requires sequence files)
+impg lace -g *.gfa -o combined.gfa --fill-gaps 2 # Fill with N's
+impg lace -g *.gfa -o combined.gfa --fill-gaps 2 --sequence-files sequence.fa # Fill with sequences
+
+# Control output compression
+impg lace -g *.gfa -o combined.gfa.gz --compress gzip
+impg lace -g *.gfa -o combined.gfa.bgz --compress bgzip
+impg lace -g *.gfa -o combined.gfa.zst --compress zstd
+
+# Use custom temporary directory
+impg lace -g *.gfa -o combined.gfa --temp-dir /tmp/lace_work
+```
+
+#### Path Name Format
+
+The command expects path names in the format:
+
+```
+NAME:START-END
+```
+
+Example: `HG002#1#chr20:1000-2000`
+
+The command uses these coordinates to:
+1. Identify which sequences belong together
+2. Order the sequences correctly
+3. Detect and handle overlaps or gaps
+
+Note: `NAME` can contain ':' characters. When parsing coordinates, the command uses the last occurrence of ':' to separate the name from the coordinate range.
+
+#### Post-processing recommendations
+
+After combining the GFA files, the resulting graph will already have compacted node IDs ranging from `1` to the total number of nodes. However, it is strongly recommended to perform post-processing steps using **[ODGI](https://github.com/pangenome/odgi)** to unchop and sort the graph.
+
+```bash
+odgi unchop -i combined.gfa -o - -t 16 | \
+    odgi sort -i - -o - -p gYs -t 16 | \
+    odgi view -i - -g > combined.final.gfa
+```
+
+If overlaps were present, and then trimmed during the merging process, it's advisable to run **[GFAffix](https://github.com/marschall-lab/GFAffix)** before the ODGI pipeline to remove redundant nodes introduced by the overlap trimming.
+
+```bash
+gfaffix combined.gfa -o combined.fix.gfa &> /dev/null
+
+odgi unchop -i combined.fix.gfa -o - -t 16 | \
+    odgi sort -i - -o - -p gYs -t 16 | \
+    odgi view -i - -g > combined.final.gfa
+```
+
 ### Index
 
 Create an IMPG index from PAF files:

README.md

@@ -53,7 +53,7 @@ impl AgcIndex {
                 if !agc.open(agc_path, true) {
                     return Err(io::Error::new(
                         io::ErrorKind::InvalidData,
-                        format!("Failed to open AGC file: {}", agc_path),
+                        format!("Failed to open AGC file: {agc_path}"),
                     ));
                 }
 
@@ -78,7 +78,7 @@ impl AgcIndex {
             for (sample, contigs) in sample_contigs {
                 for contig in contigs {
                     // Create a key that combines full contig name and sample name
-                    let key = format!("{}@{}", contig, sample);
+                    let key = format!("{contig}@{sample}");
                     index.sample_contig_to_agc.insert(key.clone(), agc_idx);
 
                     // Also insert just the full contig name if it's unique
@@ -102,7 +102,7 @@ impl AgcIndex {
                     // If short name differs from full name, also create mappings for short name
                     if short_contig != contig {
                         // Create key with short contig name and sample
-                        let short_key = format!("{}@{}", short_contig, sample);
+                        let short_key = format!("{short_contig}@{sample}");
                         index
                             .sample_contig_to_agc
                             .entry(short_key)
@@ -154,7 +154,7 @@ impl AgcIndex {
         let agc_idx = agc_idx.ok_or_else(|| {
             io::Error::new(
                 io::ErrorKind::NotFound,
-                format!("Sequence '{}' not found in any AGC file", seq_name),
+                format!("Sequence '{seq_name}' not found in any AGC file"),
             )
         })?;
 
@@ -164,8 +164,7 @@ impl AgcIndex {
             .get_contig_sequence(&sample, &contig, start, end - 1)
             .map_err(|e| {
                 io::Error::other(format!(
-                    "Failed to fetch sequence '{}@{}:{}:{}': {}",
-                    contig, sample, start, end, e
+                    "Failed to fetch sequence '{contig}@{sample}:{start}:{end}': {e}"
                 ))
             })?;
 
@@ -174,11 +173,11 @@ impl AgcIndex {
 
     pub fn get_sequence_length(&self, seq_name: &str) -> io::Result<usize> {
         let (sample, contig, agc_idx) = self.parse_query(seq_name);
-        
+
         let agc_idx = agc_idx.ok_or_else(|| {
             io::Error::new(
                 io::ErrorKind::NotFound,
-                format!("Sequence '{}' not found in any AGC file", seq_name),
+                format!("Sequence '{seq_name}' not found in any AGC file"),
             )
         })?;