orientation fix, resolves #354 (#355)

ratschlab · Nov 8, 2021 · cc68a17 · cc68a17
1 parent 1cf7a6e
commit cc68a17
Show file tree

Hide file tree

Showing 4 changed files with 47 additions and 6 deletions.
diff --git a/metagraph/src/cli/align.cpp b/metagraph/src/cli/align.cpp
@@ -386,7 +386,6 @@ int align_to_graph(Config *config) {
     for (const auto &file : files) {
         logger->trace("Align sequences from file {}", file);
         seq_io::FastaParser fasta_parser(file, config->forward_and_reverse);
-        bool is_reverse_complement = false;
 
         Timer data_reading_timer;
 
@@ -418,8 +417,7 @@ int align_to_graph(Config *config) {
                                                   config->fasta_anno_comment_delim,
                                                   true)
                             : std::string(it->name.s);
-                seq_batch.emplace_back(std::move(header), it->seq.s, is_reverse_complement);
-                is_reverse_complement ^= config->forward_and_reverse;
+                seq_batch.emplace_back(std::move(header), it->seq.s, false);
                 num_bytes_read += it->seq.l;
             }
 

diff --git a/metagraph/src/cli/config/config.cpp b/metagraph/src/cli/config/config.cpp
@@ -948,7 +948,7 @@ void Config::print_usage(const std::string &prog_name, IdentityType identity) {
             fprintf(stderr, "Usage: %s align -i <GRAPH> [options] FASTQ1 [[FASTQ2] ...]\n\n", prog_name.c_str());
 
 #if ! _PROTEIN_GRAPH
-            fprintf(stderr, "\t   --fwd-and-reverse \t\tfor each input sequence, align its reverse complement as well [off]\n");
+            fprintf(stderr, "\t   --fwd-and-reverse \t\tfor each input sequence, report a separate alignment for its reverse complement as well [off]\n");
 #endif
             fprintf(stderr, "\t   --header-comment-delim [STR]\tdelimiter for joining fasta header with comment [off]\n");
             fprintf(stderr, "\t-p --parallel [INT] \t\tuse multiple threads for computation [1]\n");

diff --git a/metagraph/src/graph/alignment/dbg_aligner.cpp b/metagraph/src/graph/alignment/dbg_aligner.cpp
@@ -79,8 +79,8 @@ ::align_batch(const std::vector<IDBGAligner::Query> &seq_batch,
             auto seeder_rc = build_seeder(reverse, !is_reverse_complement, nodes_rc);
             auto extender_rc = build_extender(reverse, aggregator);
 
-            align_both_directions(paths.get_query(false), paths.get_query(true),
-                                  *seeder, *seeder_rc, *extender, *extender_rc,
+            align_both_directions(this_query, reverse, *seeder, *seeder_rc,
+                                  *extender, *extender_rc,
                                   add_alignment, get_min_path_score);
 
             num_explored_nodes += extender_rc->num_explored_nodes();

diff --git a/metagraph/tests/graph/test_aligner.cpp b/metagraph/tests/graph/test_aligner.cpp
@@ -322,6 +322,49 @@ TYPED_TEST(DBGAlignerTest, align_straight_forward_and_reverse_complement) {
     }
 }
 
+TYPED_TEST(DBGAlignerTest, align_straight_forward_and_reverse_complement_batch) {
+    size_t k = 4;
+    std::string reference = "AGCTTCGAGGCCAA";
+    std::string query = reference;
+    reverse_complement(query.begin(), query.end());
+
+    auto graph = build_graph_batch<TypeParam>(k, { reference });
+
+    DBGAlignerConfig config(DBGAlignerConfig::dna_scoring_matrix(2, -1, -2));
+    auto config_fwd_and_rev = config;
+
+    DBGAligner<> aligner(*graph, config_fwd_and_rev);
+    std::vector<IDBGAligner::Query> query_batch { { "", query, true } };
+    aligner.align_batch(query_batch, [&](std::string_view, auto&& paths) {
+        ASSERT_EQ(1ull, paths.size());
+        auto path = paths[0];
+
+        EXPECT_EQ(query.size() - k + 1, path.size());
+        EXPECT_EQ(reference, path.get_sequence());
+        EXPECT_EQ(config.match_score(query), path.get_score());
+        EXPECT_EQ("14=", path.get_cigar().to_string());
+        EXPECT_EQ(14u, path.get_cigar().get_num_matches());
+        EXPECT_TRUE(is_exact_match(path));
+        EXPECT_EQ(0u, path.get_clipping());
+        EXPECT_EQ(0u, path.get_end_clipping());
+        EXPECT_EQ(0u, path.get_offset());
+        EXPECT_TRUE(path.is_valid(*graph, &config));
+    });
+
+    auto paths = aligner.align(query, true);
+    auto path = paths[0];
+    EXPECT_EQ(query.size() - k + 1, path.size());
+    EXPECT_EQ(reference, path.get_sequence());
+    EXPECT_EQ(config.match_score(query), path.get_score());
+    EXPECT_EQ("14=", path.get_cigar().to_string());
+    EXPECT_EQ(14u, path.get_cigar().get_num_matches());
+    EXPECT_TRUE(is_exact_match(path));
+    EXPECT_EQ(0u, path.get_clipping());
+    EXPECT_EQ(0u, path.get_end_clipping());
+    EXPECT_EQ(0u, path.get_offset());
+    EXPECT_TRUE(path.is_valid(*graph, &config));
+}
+
 #endif