seq::alg::MergePairwiseAlignments, fix gap at first position

git-svn-id: https://dng.biozentrum.unibas.ch/svn/openstructure/trunk@2567 5a81b35b-ba03-0410-adc8-b2c5c5119f08

modules/seq/alg/src/merge_pairwise_alignments.cc

@@ -50,8 +50,13 @@ void update_shifts(const AlignmentHandle& aln,
       j++;
     }
     if (shift>0) {
-      assert(i>0);
-      int res_index=s1.GetResidueIndex(i-1);
+      int res_index=0;
+      if (i == 0){
+        res_index=-1;
+      }
+      else if (i > 0) {
+        res_index=s1.GetResidueIndex(i-1);
+      }
       ShiftMap::iterator p=shifts.find(res_index);
       if (p!=shifts.end()) {
         p->second=std::max(p->second, shift);
@@ -70,13 +75,18 @@ SequenceHandle shift_reference(const SequenceHandle& ref_seq,
   String ref_str=ref_seq.GetString();
   int last=0;
   for (ShiftMap::const_iterator i=shifts.begin(),e=shifts.end();i!=e; ++i) {
-    new_sequence << ref_str.substr(last, i->first-last+1) 
+    if(i->first == -1){
+      new_sequence << String(i->second,'-');
+    }
+    else{
+      new_sequence << ref_str.substr(last, i->first-last+1)
                  << String(i->second, '-');
-    last=i->first+1;
+      last=i->first+1;
+    }
   }
   new_sequence << ref_str.substr(last);
   SequenceHandle s=CreateSequence(ref_seq.GetName(), 
-                                            new_sequence.str());                                               
+                                            new_sequence.str());
   if (ref_seq.HasAttachedView())
     s.AttachView(ref_seq.GetAttachedView());
   s.SetSequenceOffset(s.GetSequenceOffset());
@@ -96,15 +106,19 @@ SequenceHandle realign_sequence(const AlignmentHandle& aln,
       shift++;
       ++j;
     }
-    if (i>0 && s1.GetOneLetterCode(i-1)!='-') {
-      ShiftMap::const_iterator p=shifts.find(s1.GetResidueIndex(i-1));
-      if (p!=shifts.end()) {
-        int d=p->second-shift;
-        assert(d>=0);
-        new_sequence << String(d, '-');
-      }      
+    ShiftMap::const_iterator p=shifts.end();
+    if (i == 0){
+      p=shifts.find(-1);
+    }
+    else if (i>0 && s1.GetOneLetterCode(i-1)!='-') {
+      p=shifts.find(s1.GetResidueIndex(i-1));
+    }
+    if (p!=shifts.end()) {
+      int d=p->second-shift;
+      assert(d>=0);
+      new_sequence << String(d, '-');
     }
-    new_sequence << s2.GetOneLetterCode(i);    
+    new_sequence << s2.GetOneLetterCode(i);
   }
   SequenceHandle s=CreateSequence(s2.GetName(), new_sequence.str());
   if (s2.HasAttachedView())

modules/seq/alg/tests/test_merge_pairwise_alignments.cc

@@ -60,8 +60,8 @@ BOOST_AUTO_TEST_CASE(merge_pairwise_alignments_one)
 BOOST_AUTO_TEST_CASE(merge_pairwise_alignments_two) 
 {
   SequenceHandle ref=CreateSequence("REF", "abcdefghijklmn");
-  SequenceHandle s1=CreateSequence("S1", "abcd---efghijklmn");
-  SequenceHandle s2=CreateSequence("S2", "abcdxyzefghijklmn");
+  SequenceHandle s1=CreateSequence("S1",   "abcd---efghijklmn");
+  SequenceHandle s2=CreateSequence("S2",   "abcdxyzefghijklmn");
   AlignmentHandle aln1=CreateAlignment();
   aln1.AddSequence(s1);
   aln1.AddSequence(s2);
@@ -82,4 +82,29 @@ BOOST_AUTO_TEST_CASE(merge_pairwise_alignments_two)
   BOOST_CHECK_EQUAL(seqs[2].GetString(), "abcd--yefghijklmn");
 }
 
+
+BOOST_AUTO_TEST_CASE(merge_pairwise_alignments_three)
+{
+  SequenceHandle ref=CreateSequence("REF", "abcdefghijk");
+  SequenceHandle s1=CreateSequence("S1",   "--abcdefghijk");
+  SequenceHandle s2=CreateSequence("S2",   "xyabcdefghijk");
+  AlignmentHandle aln1=CreateAlignment();
+  aln1.AddSequence(s1);
+  aln1.AddSequence(s2);
+
+  SequenceHandle s3=CreateSequence("S1", "-abcdefghijk");
+  SequenceHandle s4=CreateSequence("S2", "zabcdefghijk");
+
+  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[1].GetString(), "xyabcdefghijk");
+  BOOST_CHECK_EQUAL(seqs[2].GetString(), "-zabcdefghijk");
+}
 BOOST_AUTO_TEST_SUITE_END()