Support outputting alignments in GAF format

metagraph/src/cli/align.cpp

@@ -287,18 +287,7 @@ std::string format_alignment(std::string_view header,
                              const DeBruijnGraph &graph,
                              const Config &config) {
     std::string sout;
-    if (!config.output_json) {
-        sout += fmt::format("{}\t{}", header, paths.get_query());
-        if (paths.empty()) {
-            sout += fmt::format("\t*\t*\t{}\t*\t*\t*", config.alignment_min_path_score);
-        } else {
-            for (const auto &path : paths.data()) {
-                sout += fmt::format("\t{}", path);
-            }
-        }
-
-        sout += "\n";
-    } else {
+    if (config.output_json) {
         Json::StreamWriterBuilder builder;
         builder["indentation"] = "";
 
@@ -319,6 +308,58 @@ std::string format_alignment(std::string_view header,
 
             sout += fmt::format("{}\n", Json::writeString(builder, json_line));
         }
+    } else if (config.output_gaf) {
+        std::string base = fmt::format("{}\t{}\t", header, paths.get_query().size());
+        for (size_t i = 0; i < paths.size(); ++i) {
+            const auto &path = paths[i];
+            const auto &cigar = path.get_cigar();
+            const auto &nodes = path.get_nodes();
+            sout += base;
+            if (path.empty()) {
+                sout += "*\t*\t*\t*\t*\t*\t*\t*\t*\t*\n";
+                continue;
+            }
+
+            size_t begin = 0;
+            size_t end = 0;
+
+            if (path.get_orientation()) {
+                begin = cigar.get_end_clipping();
+                end = paths.get_query().size() - cigar.get_clipping();
+            } else {
+                begin = cigar.get_clipping();
+                end = paths.get_query().size() - cigar.get_end_clipping();
+            }
+
+            sout += fmt::format(
+                "{}\t{}\t+\t{}\t{}\t{}\t{}\t{}\t{}\t255\tAS:i:{}\tcg:Z:{}\tMD:Z:{}",
+                begin, end,
+                path.get_orientation()
+                    ? fmt::format("<{}", fmt::join(nodes.rbegin(), nodes.rend(), ">"))
+                    : fmt::format(">{}", fmt::join(nodes, ">")),
+                path.get_sequence().size() + path.get_offset(),
+                path.get_orientation() ? path.get_offset() : 0,
+                !path.get_orientation() ? path.get_offset() : 0,
+                cigar.get_num_matches(),
+                path.get_sequence().size(),
+                path.get_score(),
+                cigar.to_string(),
+                cigar.to_md_string(path.get_sequence())
+            );
+
+            sout += "\n";
+        }
+    } else {
+        sout += fmt::format("{}\t{}", header, paths.get_query());
+        if (paths.empty()) {
+            sout += fmt::format("\t*\t*\t{}\t*\t*\t*", config.alignment_min_path_score);
+        } else {
+            for (const auto &path : paths.data()) {
+                sout += fmt::format("\t{}", path);
+            }
+        }
+
+        sout += "\n";
     }
 
     return sout;

metagraph/src/cli/config/config.cpp

@@ -338,6 +338,8 @@ Config::Config(int argc, char *argv[]) {
             output_compacted = true;
         } else if (!strcmp(argv[i], "--json")) {
             output_json = true;
+        } else if (!strcmp(argv[i], "--gaf")) {
+            output_gaf = true;
         } else if (!strcmp(argv[i], "--unitigs")) {
             to_fasta = true;
             unitigs = true;
@@ -1005,6 +1007,7 @@ if (advanced) {
             fprintf(stderr, "Available options for alignment:\n");
             fprintf(stderr, "\t-o --outfile-base [STR]\t\t\t\tbasename of output file []\n");
             fprintf(stderr, "\t   --json \t\t\t\t\toutput alignment in JSON format [off]\n");
+            fprintf(stderr, "\t   --gaf \t\t\t\t\toutput alignment in GAF format [off]\n");
 if (advanced) {
             fprintf(stderr, "\t   --align-only-forwards \t\t\tdo not align backwards from a seed on basic-mode graphs [off]\n");
 }

metagraph/src/cli/config/config.hpp

@@ -61,6 +61,7 @@ class Config {
     bool align_only_forwards = false;
     bool filter_by_kmer = false;
     bool output_json = false;
+    bool output_gaf = false;
     bool aggregate_columns = false;
     bool coordinates = false;
     bool advanced = false;

metagraph/src/graph/alignment/aligner_cigar.cpp

@@ -93,6 +93,46 @@ std::string Cigar::to_string() const {
     return cigar_string;
 }
 
+std::string Cigar::to_md_string(std::string_view reference) const {
+    std::string md_string;
+    auto ref_it = reference.begin();
+
+    size_t match_count = 0;
+    for (const auto &[op, num] : cigar_) {
+        switch (op) {
+            case CLIPPED:
+            case INSERTION:
+            case NODE_INSERTION: {} break;
+            case MATCH: {
+                match_count += num;
+                ref_it += num;
+            } break;
+            case MISMATCH: {
+                if (match_count) {
+                    md_string += std::to_string(match_count);
+                    match_count = 0;
+                }
+                md_string += std::string(ref_it, ref_it + num);
+                ref_it += num;
+            } break;
+            case DELETION: {
+                if (match_count) {
+                    md_string += std::to_string(match_count);
+                    match_count = 0;
+                }
+                md_string += "^" + std::string(ref_it, ref_it + num);
+                ref_it += num;
+            }
+        }
+    }
+
+    if (match_count)
+        md_string += std::to_string(match_count);
+
+    assert(ref_it == reference.end());
+    return md_string;
+}
+
 void Cigar::append(Operator op, LengthType num) {
     if (!num)
         return;

metagraph/src/graph/alignment/aligner_cigar.hpp

@@ -40,6 +40,7 @@ class Cigar {
     bool empty() const { return cigar_.empty(); }
 
     std::string to_string() const;
+    std::string to_md_string(std::string_view reference) const;
 
     void append(Operator op, LengthType num = 1);
     void append(Cigar&& other);