diff --git a/modules/seq/alg/src/merge_pairwise_alignments.cc b/modules/seq/alg/src/merge_pairwise_alignments.cc index f7a7c000ecebf9cc96711b59dff976a6c9620816..664f6ccc6315da166fd169daa20951d4f1b32367 100644 --- a/modules/seq/alg/src/merge_pairwise_alignments.cc +++ b/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()) diff --git a/modules/seq/alg/tests/test_merge_pairwise_alignments.cc b/modules/seq/alg/tests/test_merge_pairwise_alignments.cc index 37d58e7073eebf2726fd5b8e170bb9096fec78bf..4e0e4d67a7bccc504df611dad5c1bad826c38cab 100644 --- a/modules/seq/alg/tests/test_merge_pairwise_alignments.cc +++ b/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()