More fixes, docs and unitests for MergePairwiseAlignment

modules/seq/alg/doc/seqalg.rst

@@ -22,6 +22,10 @@
    - The first sequence of each pairwise alignments corresponds to the reference 
      sequence. Apart from the presence of gaps, these two sequences must be 
      completely identical.
+
+   - If the reference sequence has an offset, the first sequence of each pairwise alignment 
+     must have the same offset. This offset is inherited by the first sequence of the final
+     output alignment.
    
    - The resulting multiple sequence alignment is by no means optimal. For 
      better results, consider using a multiple-sequence alignment program such 

modules/seq/alg/src/merge_pairwise_alignments.cc

@@ -90,7 +90,7 @@ SequenceHandle shift_reference(const ConstSequenceHandle& ref_seq,
                                   new_sequence.str());
   if (ref_seq.HasAttachedView())
     s.AttachView(ref_seq.GetAttachedView());
-  s.SetOffset(s.GetOffset());
+  s.SetOffset(ref_seq.GetOffset());
   return s;
 }
 
@@ -134,7 +134,23 @@ SequenceHandle realign_sequence(const AlignmentHandle& aln,
 AlignmentHandle MergePairwiseAlignments(const AlignmentList& pairwise_alns,
                                         const ConstSequenceHandle& ref_seq)
 {
-
+   
+  int ref_offset=ref_seq.GetOffset();
+  int alignment_counter=0;
+  for (AlignmentList::const_iterator i=pairwise_alns.begin(), e=pairwise_alns.end(); i!=e; ++i) {
+    AlignmentHandle ali_handle=*i;
+    ConstSequenceHandle seq_handle=ali_handle.GetSequence(0);
+    if (seq_handle.GetOffset()!=ref_offset) {
+      String error_string;
+      std::ostringstream error_ss(error_string);
+      error_ss << "The offset of the first sequence in alignment " << alignment_counter << 
+        "is not identical to the offset of the reference sequence." << std::endl;
+      throw IntegrityError(error_string);
+    }
+    alignment_counter+=1;
+  }
+  
+  
   ShiftMap shifts;
   for (AlignmentList::const_iterator i=pairwise_alns.begin(),
        e=pairwise_alns.end(); i!=e; ++i) {

modules/seq/alg/tests/test_merge_pairwise_alignments.cc

@@ -24,7 +24,7 @@
 #include <boost/test/unit_test.hpp>
 
 #include <ost/seq/alg/merge_pairwise_alignments.hh>
-
+#include <ost/integrity_error.hh>
 
 
 using namespace ost;
@@ -133,4 +133,60 @@ BOOST_AUTO_TEST_CASE(merge_pairwise_alignments_four)
   BOOST_CHECK_EQUAL(seqs[2].GetString(), "zabcdefghijk--");
 }
 
+BOOST_AUTO_TEST_CASE(merge_pairwise_alignments_five)
+{
+  SequenceHandle ref=CreateSequence("REF", "abcdefghijk");
+  ref.SetOffset(4);
+  SequenceHandle s1=CreateSequence("S1",   "abcdefghijk--");
+  SequenceHandle s2=CreateSequence("S2",   "abcdefghijkxy");
+  s1.SetOffset(4);
+  AlignmentHandle aln1=CreateAlignment();
+  aln1.AddSequence(s1); 
+  aln1.AddSequence(s2);
+
+  SequenceHandle s3=CreateSequence("S1", "-abcdefghijk");
+  SequenceHandle s4=CreateSequence("S2", "zabcdefghijk");
+  s3.SetOffset(4);
+  
+  AlignmentHandle aln2=CreateAlignment();
+  aln2.AddSequence(s3);
+  aln2.AddSequence(s4);
+  AlignmentList l;
+  l.push_back(aln1);
+  l.push_back(aln2);
+  AlignmentHandle m=alg::MergePairwiseAlignments(l, ref);
+  ConstSequenceList seqs=m.GetSequences();  
+  
+  BOOST_CHECK_EQUAL(seqs[0].GetString(), "-abcdefghijk--");
+  BOOST_CHECK_EQUAL(seqs[0].GetOffset(), 4);
+  BOOST_CHECK_EQUAL(seqs[1].GetString(), "-abcdefghijkxy");
+  BOOST_CHECK_EQUAL(seqs[2].GetString(), "zabcdefghijk--");
+}
+
+BOOST_AUTO_TEST_CASE(merge_pairwise_alignments_six)
+{
+  SequenceHandle ref=CreateSequence("REF", "abcdefghijk");
+  ref.SetOffset(4);
+  SequenceHandle s1=CreateSequence("S1",   "abcdefghijk--");
+  SequenceHandle s2=CreateSequence("S2",   "abcdefghijkxy");
+  s1.SetOffset(5);
+  AlignmentHandle aln1=CreateAlignment();
+  aln1.AddSequence(s1); 
+  aln1.AddSequence(s2);
+
+  SequenceHandle s3=CreateSequence("S1", "-abcdefghijk");
+  SequenceHandle s4=CreateSequence("S2", "zabcdefghijk");
+  s3.SetOffset(4);
+  
+  AlignmentHandle aln2=CreateAlignment();
+  aln2.AddSequence(s3);
+  aln2.AddSequence(s4);
+  AlignmentList l;
+  l.push_back(aln1);
+  l.push_back(aln2);
+  BOOST_CHECK_THROW (alg::MergePairwiseAlignments(l, ref),IntegrityError);
+  
+}
+
+
 BOOST_AUTO_TEST_SUITE_END()