diff --git a/modules/seq/alg/doc/seqalg.rst b/modules/seq/alg/doc/seqalg.rst index cf3587d06c21d3041cbe350bad899dde7dfdef82..79e6ca62bb3d97082b37fcbb9afee719e12b276f 100644 --- a/modules/seq/alg/doc/seqalg.rst +++ b/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 diff --git a/modules/seq/alg/src/merge_pairwise_alignments.cc b/modules/seq/alg/src/merge_pairwise_alignments.cc index d5da17127fc3e5233552a721d12acdaa1f1eb0f5..9d60859a07b99b7dda1724f54306f91a7605690b 100644 --- a/modules/seq/alg/src/merge_pairwise_alignments.cc +++ b/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) { diff --git a/modules/seq/alg/tests/test_merge_pairwise_alignments.cc b/modules/seq/alg/tests/test_merge_pairwise_alignments.cc index f89715adeccc627198420c1fc80826502db617a6..851f5c9f163911f1166492a25fd94ebcf4c5ea5a 100644 --- a/modules/seq/alg/tests/test_merge_pairwise_alignments.cc +++ b/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()