Fix support to GFA2 & also support GFA1, fixes #59

Author: Guilucand

crates/api/example/src/main.rs

@@ -27,7 +27,6 @@ fn main() {
                 panic!("Unrecoverable error: {}", msg);
             }
         }),
-        gfa_output: false,
         disk_optimization_level: 5,
     })
     .unwrap();
@@ -63,7 +62,7 @@ fn main() {
             true,
             1,
             ExtraElaboration::UnitigLinks,
-            false,
+            None,
             5,
         )
         .unwrap();

crates/api/src/lib.rs

@@ -9,7 +9,7 @@ use config::{MAX_BUCKET_CHUNK_SIZE, MIN_BUCKET_CHUNK_SIZE};
 pub use ggcat_logging::MessageLevel;
 use ggcat_logging::UnrecoverableErrorLogging;
 use io::concurrent::structured_sequences::fasta::FastaWriterWrapper;
-use io::concurrent::structured_sequences::gfa::GFAWriterWrapper;
+use io::concurrent::structured_sequences::gfa::{GFAWriterWrapperV1, GFAWriterWrapperV2};
 use io::concurrent::structured_sequences::StructuredSequenceBackendWrapper;
 use io::sequences_stream::fasta::FastaFileSequencesStream;
 use io::sequences_stream::GenericSequencesStream;
@@ -59,11 +59,18 @@ pub mod debug {
     pub static BUCKETS_COUNT_LOG_FORCE: Mutex<Option<usize>> = Mutex::new(None);
 }
 
+#[derive(Clone, Copy, Debug)]
 pub enum LoggingMode {
     Log,
     ForceStdout,
 }
 
+#[derive(Clone, Copy, Debug)]
+pub enum GfaVersion {
+    V1,
+    V2,
+}
+
 /// Main config of GGCAT. This config is global and should be passed to GGCATInstance::create
 pub struct GGCATConfig {
     /// Directory for temporary files
@@ -90,9 +97,6 @@ pub struct GGCATConfig {
     /// The messages callback, if present, no output will be automatically written to stdout
     pub messages_callback: Option<fn(MessageLevel, &str)>,
 
-    /// Output GFA format instead of FASTA
-    pub gfa_output: bool,
-
     /// Sets the level of disk usage reduction optimization
     pub disk_optimization_level: u32,
 }
@@ -231,7 +235,7 @@ impl GGCATInstance {
 
         extra_elab: ExtraElaboration,
 
-        gfa_output: bool,
+        gfa_output_version: Option<GfaVersion>,
 
         disk_optimization_level: u32,
     ) -> anyhow::Result<PathBuf> {
@@ -249,14 +253,14 @@ impl GGCATInstance {
             NonColoredManager::dynamic_dispatch_id()
         };
 
-        if gfa_output && colors {
+        if gfa_output_version.is_some() && colors {
             anyhow::bail!("GFA output is not supported with colors");
         }
 
-        let output_mode = if gfa_output {
-            GFAWriterWrapper::dynamic_dispatch_id()
-        } else {
-            FastaWriterWrapper::dynamic_dispatch_id()
+        let output_mode = match gfa_output_version {
+            None => FastaWriterWrapper::dynamic_dispatch_id(),
+            Some(GfaVersion::V1) => GFAWriterWrapperV1::dynamic_dispatch_id(),
+            Some(GfaVersion::V2) => GFAWriterWrapperV2::dynamic_dispatch_id(),
         };
 
         let temp_dir = create_tempdir(self.0.temp_dir.clone());

crates/assembler/src/lib.rs

@@ -19,7 +19,8 @@ use config::{
 use hashes::HashFunctionFactory;
 use io::concurrent::structured_sequences::binary::StructSeqBinaryWriter;
 use io::concurrent::structured_sequences::fasta::FastaWriterWrapper;
-use io::concurrent::structured_sequences::gfa::GFAWriterWrapper;
+use io::concurrent::structured_sequences::gfa::GFAWriterWrapperV1;
+use io::concurrent::structured_sequences::gfa::GFAWriterWrapperV2;
 use io::concurrent::structured_sequences::{
     IdentSequenceWriter, StructuredSequenceBackend, StructuredSequenceBackendInit,
     StructuredSequenceBackendWrapper, StructuredSequenceWriter,
@@ -95,7 +96,8 @@ fn get_writer<
     colors::non_colored::NonColoredManager,
 ], OutputMode = [
     FastaWriterWrapper,
-    #[cfg(not(feature = "devel-build"))] GFAWriterWrapper
+    #[cfg(not(feature = "devel-build"))] GFAWriterWrapperV1,
+    #[cfg(not(feature = "devel-build"))] GFAWriterWrapperV2
 ])]
 pub fn run_assembler<
     MergingHash: HashFunctionFactory,

crates/assembler/src/pipeline/maximal_unitig_links.rs

@@ -172,6 +172,7 @@ pub fn build_maximal_unitigs_links<
                                         index,
                                         MaximalUnitigPosition::Beginning,
                                         first_hash.is_forward(),
+                                        0,
                                     ),
                                 );
 
@@ -188,6 +189,7 @@ pub fn build_maximal_unitigs_links<
                                             index,
                                             MaximalUnitigPosition::Beginning,
                                             !first_hash.is_forward(),
+                                            0,
                                         ),
                                     );
                                 }
@@ -204,6 +206,7 @@ pub fn build_maximal_unitigs_links<
                                         index,
                                         MaximalUnitigPosition::Ending,
                                         !last_hash.is_forward(),
+                                        read_len as u64 - k as u64 + 1,
                                     ),
                                 );
 
@@ -220,6 +223,7 @@ pub fn build_maximal_unitigs_links<
                                             index,
                                             MaximalUnitigPosition::Ending,
                                             last_hash.is_forward(),
+                                            read_len as u64 - k as u64 + 1,
                                         ),
                                     );
                                 }
@@ -310,6 +314,7 @@ pub fn build_maximal_unitigs_links<
                             .map(|v| {
                                 MaximalUnitigIndex::new(
                                     v.entry(),
+                                    v.overlap_start(),
                                     MaximalUnitigFlags::new_direction(
                                         val.position() == MaximalUnitigPosition::Beginning,
                                         v.position() == MaximalUnitigPosition::Ending,

crates/assembler/src/pipeline/maximal_unitig_links/maximal_hash_entry.rs

@@ -19,6 +19,7 @@ pub enum MaximalUnitigPosition {
 pub struct MaximalHashEntry<H: Copy> {
     pub hash: H,
     encoded: u64,
+    overlap_start: u64,
 }
 
 impl<H: Copy> MaximalHashEntry<H> {
@@ -31,6 +32,7 @@ impl<H: Copy> MaximalHashEntry<H> {
         entry: u64,
         position: MaximalUnitigPosition,
         direction_forward: bool,
+        overlap_start: u64,
     ) -> Self {
         Self {
             hash,
@@ -40,13 +42,18 @@ impl<H: Copy> MaximalHashEntry<H> {
                     MaximalUnitigPosition::Beginning => 0,
                 }) << Self::POSITION_OFFSET)
                 | ((if direction_forward { 1 } else { 0 }) << Self::DIRECTION_OFFSET),
+            overlap_start,
         }
     }
 
     pub fn entry(&self) -> u64 {
         self.encoded >> Self::ENTRY_OFFSET
     }
 
+    pub fn overlap_start(&self) -> u64 {
+        self.overlap_start
+    }
+
     pub fn position(&self) -> MaximalUnitigPosition {
         if (self.encoded >> Self::POSITION_OFFSET) & 0x1 == 0 {
             MaximalUnitigPosition::Beginning

crates/assembler/src/pipeline/maximal_unitig_links/maximal_unitig_index.rs

@@ -44,6 +44,8 @@ impl MaximalUnitigFlags {
 #[derive(Copy, Clone, Eq)]
 pub struct MaximalUnitigIndex {
     index: u64,
+    // The non reverse-complemented start of the link overlap (used for gfa2)
+    overlap_start: u64,
     pub flags: MaximalUnitigFlags,
 }
 
@@ -63,12 +65,18 @@ impl HasEmptyExtraBuffer for MaximalUnitigIndex {}
 impl SequenceExtraData for MaximalUnitigIndex {
     fn decode_extended(_: &mut (), reader: &mut impl Read) -> Option<Self> {
         let index = decode_varint(|| reader.read_u8().ok())?;
+        let overlap_start = decode_varint(|| reader.read_u8().ok())?;
         let flags = reader.read_u8().ok()?;
-        Some(MaximalUnitigIndex::new(index, MaximalUnitigFlags(flags)))
+        Some(MaximalUnitigIndex::new(
+            index,
+            overlap_start,
+            MaximalUnitigFlags(flags),
+        ))
     }
 
     fn encode_extended(&self, _: &(), writer: &mut impl Write) {
         encode_varint(|b| writer.write_all(b).ok(), self.index() as u64).unwrap();
+        encode_varint(|b| writer.write_all(b).ok(), self.overlap_start).unwrap();
         writer.write_all(&[self.flags.0]).unwrap();
     }
 
@@ -88,8 +96,12 @@ impl Debug for MaximalUnitigIndex {
 
 impl MaximalUnitigIndex {
     #[inline]
-    pub fn new(index: u64, flags: MaximalUnitigFlags) -> Self {
-        Self { index, flags }
+    pub fn new(index: u64, overlap_start: u64, flags: MaximalUnitigFlags) -> Self {
+        Self {
+            index,
+            overlap_start,
+            flags,
+        }
     }
 
     #[inline]
@@ -148,6 +160,7 @@ impl BucketItemSerializer for MaximalUnitigLinkSerializer {
 
         for entry in entries {
             encode_varint(|b| bucket.write_all(b), entry.index() as u64).unwrap();
+            encode_varint(|b| bucket.write_all(b), entry.overlap_start).unwrap();
             bucket.push(entry.flags.0);
         }
     }
@@ -165,8 +178,13 @@ impl BucketItemSerializer for MaximalUnitigLinkSerializer {
         let start = read_buffer.len();
         for _i in 0..len {
             let index = decode_varint(|| stream.read_u8().ok())?;
+            let overlap_start = decode_varint(|| stream.read_u8().ok())?;
             let flags = stream.read_u8().ok()?;
-            read_buffer.push(MaximalUnitigIndex::new(index, MaximalUnitigFlags(flags)));
+            read_buffer.push(MaximalUnitigIndex::new(
+                index,
+                overlap_start,
+                MaximalUnitigFlags(flags),
+            ));
         }
 
         Some(MaximalUnitigLink::new(entry, VecSlice::new(start, len)))
@@ -272,26 +290,61 @@ impl IdentSequenceWriter for DoubleMaximalUnitigLinks {
     }
 
     #[allow(unused_variables)]
-    fn write_as_gfa(
+    fn write_as_gfa<const VERSION: u32>(
         &self,
         k: u64,
         index: u64,
+        length: u64,
         stream: &mut impl Write,
         extra_buffer: &Self::TempBuffer,
     ) {
         for entries in &self.links {
             let entries = entries.entries.get_slice(extra_buffer);
-            for entry in entries {
-                writeln!(
-                    stream,
-                    "L\t{}\t{}\t{}\t{}\t{}M",
-                    index,
-                    if entry.flags.flip_current() { "-" } else { "+" },
-                    entry.index,
-                    if entry.flags.flip_other() { "-" } else { "+" },
-                    k - 1
-                )
-                .unwrap();
+            if VERSION == 1 {
+                for entry in entries {
+                    // L <index> < +/- > <other_index> < +/- > <overlap>
+                    writeln!(
+                        stream,
+                        "L\t{}\t{}\t{}\t{}\t{}M",
+                        index,
+                        if entry.flags.flip_current() { '-' } else { '+' },
+                        entry.index,
+                        if entry.flags.flip_other() { '-' } else { '+' },
+                        k - 1
+                    )
+                    .unwrap();
+                }
+            } else if VERSION == 2 {
+                for entry in entries {
+                    let (b1, e1, is_end1) = if entry.flags.flip_current() {
+                        (0, k - 1, false)
+                    } else {
+                        (length - k + 1, length, true)
+                    };
+
+                    let (b2, e2, is_end2) = (
+                        entry.overlap_start,
+                        entry.overlap_start + k - 1,
+                        entry.overlap_start > 0,
+                    );
+
+                    // E * <index>< +/- > <other_index>< +/- > <b1> <e1> <b2> <e2>
+                    writeln!(
+                        stream,
+                        "E\t*\t{}{}\t{}{}\t{}\t{}{}\t{}\t{}{}",
+                        index,
+                        if entry.flags.flip_current() { '-' } else { '+' },
+                        entry.index,
+                        if entry.flags.flip_other() { '-' } else { '+' },
+                        b1,
+                        e1,
+                        if is_end1 { "$" } else { "" },
+                        b2,
+                        e2,
+                        if is_end2 { "$" } else { "" }
+                    )
+                    .unwrap();
+                }
             }
         }
     }

crates/capi/build.rs

@@ -11,7 +11,7 @@ fn main() {
             std::env::var("OUT_DIR").unwrap(),
             "/cxxbridge/include/ggcat-cpp-bindings/src/lib.rs.h"
         ),
-        "ggcat-cpp-api/src/ggcat-cpp-bindings.hh",
+        "ggcat-cpp-api/include/ggcat-cpp-bindings.hh",
     )
     .unwrap();