From 9c9d1ed85864b35b2d88f35b4bf5c963f172cc14 Mon Sep 17 00:00:00 2001 From: Harun Mustafa Date: Mon, 7 Feb 2022 15:44:29 +0100 Subject: [PATCH] Support outputting alignments in GAF format --- metagraph/src/cli/align.cpp | 65 +++++++++++++++---- metagraph/src/cli/config/config.cpp | 3 + metagraph/src/cli/config/config.hpp | 1 + .../src/graph/alignment/aligner_cigar.cpp | 40 ++++++++++++ .../src/graph/alignment/aligner_cigar.hpp | 1 + 5 files changed, 98 insertions(+), 12 deletions(-) diff --git a/metagraph/src/cli/align.cpp b/metagraph/src/cli/align.cpp index 2f7038c2da..3e5ca8b8d0 100644 --- a/metagraph/src/cli/align.cpp +++ b/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; diff --git a/metagraph/src/cli/config/config.cpp b/metagraph/src/cli/config/config.cpp index 44b0cfd894..f6701bc089 100644 --- a/metagraph/src/cli/config/config.cpp +++ b/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"); } diff --git a/metagraph/src/cli/config/config.hpp b/metagraph/src/cli/config/config.hpp index 14693a99f9..02f4e69b2d 100644 --- a/metagraph/src/cli/config/config.hpp +++ b/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; diff --git a/metagraph/src/graph/alignment/aligner_cigar.cpp b/metagraph/src/graph/alignment/aligner_cigar.cpp index 91373ed61c..7afc68f245 100644 --- a/metagraph/src/graph/alignment/aligner_cigar.cpp +++ b/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; diff --git a/metagraph/src/graph/alignment/aligner_cigar.hpp b/metagraph/src/graph/alignment/aligner_cigar.hpp index 2ec792cf10..73a846fd1e 100644 --- a/metagraph/src/graph/alignment/aligner_cigar.hpp +++ b/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);