From b4e92f73b6b4aafaaa7ac58f028cda980f8bd8ce Mon Sep 17 00:00:00 2001
From: stefan <stefan@5a81b35b-ba03-0410-adc8-b2c5c5119f08>
Date: Thu, 15 Jul 2010 14:07:25 +0000
Subject: [PATCH] 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
---
.../seq/alg/src/merge_pairwise_alignments.cc | 40 +++++++++++++------
.../tests/test_merge_pairwise_alignments.cc | 29 +++++++++++++-
2 files changed, 54 insertions(+), 15 deletions(-)
diff --git a/modules/seq/alg/src/merge_pairwise_alignments.cc b/modules/seq/alg/src/merge_pairwise_alignments.cc
index f7a7c000e..664f6ccc6 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 37d58e707..4e0e4d67a 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()
--
GitLab