Handle subseq when the parent is on the minus strand

diff --git a/alg/seq_span.cpp b/alg/seq_span.cpp
index bf9a3dc3664bb72d40a613a28521795db1265402..1eb5e9344982cce4cf04b1b6ce364a14149de3a9 100644 (file)
--- a/alg/seq_span.cpp
+++ b/alg/seq_span.cpp
@@ -113,6 +113,14 @@ SeqSpan::SeqSpan(const SeqSpanRef parent_,
     throw sequence_invalid_strand("unrecognized strand identifier");
     break;
   }
+  
+  // Ack the complexity increases!
+  // If our parent is on the minus strand, we need to adjust the start
+  // and count to look like we're selecting from the right side of the 
+  // parent sequence (AKA the start of the minus strand)
+  if (parent and parent->strand() == MinusStrand) {
+    seq_start = parent->start() + parent->size() - (start_ + seq_count);
+  }
 }
 
 //////

diff --git a/alg/test/test_seq_span.cpp b/alg/test/test_seq_span.cpp
index 4989b6cbc8dc37159ed8604da6b71990a8b1e199..b646c1fe41480573598014f59476ba04731338dc 100644 (file)
--- a/alg/test/test_seq_span.cpp
+++ b/alg/test/test_seq_span.cpp
@@ -288,6 +288,28 @@ BOOST_AUTO_TEST_CASE( seqspan_global_mutable_start_stop )
   BOOST_CHECK_EQUAL( s3->size(), 1);
 }
 
+BOOST_AUTO_TEST_CASE( seqspan_global_mutable_start_stop_minus_strand )
+{
+  std::string seq_string("AAAAGCTA");
+  SeqSpanRef s1(new SeqSpan(seq_string));
+
+  SeqSpanRef s2 = s1->subseq(2,3, SeqSpan::MinusStrand);
+  BOOST_CHECK_EQUAL( s2->start(), 2);
+  BOOST_CHECK_EQUAL( s2->stop(), 2+3);
+  BOOST_CHECK_EQUAL( s2->size(), 3);
+  BOOST_CHECK_EQUAL( s2->sequence(), "CTT");
+  
+  SeqSpanRef s3 = s2->subseq(1,2, SeqSpan::SameStrand);
+  BOOST_CHECK_EQUAL(s3->sequence(), "TT");
+  
+  // Could also argue that it should be CT
+  // if you assume that the locations are all relative to the global sequence
+  // and are then reverse complemented
+  
+  s2->setStart(1);
+  BOOST_CHECK_EQUAL( s2->sequence(), "CTTT");
+}
+
 BOOST_AUTO_TEST_CASE( seqspan_parent_mutable_start_stop )
 {
   std::string seq_string("AAGGCCTT");