banded alignment gives spurious gaps for 1bp sequences

[mchaisso]

This isn't a common case, but it's worth bringing up as a base-case issue. Passing a single nucleotide to banded global alignment can give spurious alignments.

#include "seqan/align.h"
#include <string>
#include <iostream>
using namespace std;
typedef seqan::String<seqan::Dna> TSequence;                 // sequence type
typedef seqan::Align<TSequence, seqan::ArrayGaps> TSeqAlign;   
typedef seqan::Row<TSeqAlign>::Type TSeqRow;                 // gapped sequence typ



int main(int argc, char* argv[]) {
	string read="T";
	string ref="T";
	TSeqAlign align;
	seqan::resize(seqan::rows(align), 2);
	seqan::assignSource(seqan::row(align, 0), read.c_str());
	seqan::assignSource(seqan::row(align, 1), ref.c_str());
	int k=7;
	int score=0;
// affine banded alignment
	score = seqan::globalAlignment(align, seqan::Score<int, seqan::Simple>(4, -4, -20, 0), -k, k);
	cout << score<< endl;
	cout << align << endl;
// standard affine
	score = seqan::globalAlignment(align, seqan::Score<int, seqan::Simple>(4, -4, -20, 0));
	cout << score<< endl;
	cout << align << endl;
// standard banded
	score = seqan::globalAlignment(align, seqan::Score<int, seqan::Simple>(4, -4, -4), -k, k);
	cout << score<< endl;
	cout << align << endl;
	return 0;
}

Output:
1073741824
      0    
        T-
          
        -T



4
      0   
        T
        |
        T



-1073741824
      0    
        T-
          
        -T

[marehr]

Hi @mchaisso,

sorry for coming back to you so late, but it is currently vacation time and our group kinda focuses on version 3 of SeqAn.

You seem to have found a neat little edge case. I imagine that the 2bp case works correctly?
Could you try to use seqan::String<char> directly for read and ref without the whole .c_str() conversion?

Could you clarify what you mean by spurious?

Do you mean that affine, banded, affine-banded behave inconsistently (and wrong)? Or do you mean that multiple executions of the same binary produces different output?

Anyway, I put @rrahn which is the maintainer of the alignment code and should fairly easy find the issue. (Who is on vacation 😄 )

[mchaisso]

Hi Marcel, No hurry on this- I've actually switched to my own k-band alignment (but would be happy to switch back if there is SSE support). I believe the 2 bps case works, but I have to confirm. The spurious alignment is because the optimal alignment of align("T", "T") should be one match, but it is an insertion and a deletion. All methods are deterministic (no change in output). The affine behaves correctly, but all invocations of banded (banded and affine-banded) are incorrect. Best, -Mark

…

On Tue, Jul 17, 2018 at 2:33 AM Marcel ***@***.***> wrote: Hi @mchaisso <https://github.com/mchaisso>, sorry for coming back to you so late, but it is currently vacation time and our group kinda focuses on version 3 of SeqAn. You seem to have found a neat little edge case. I imagine that the 2bp case works correctly? Could you try to use seqan::String<char> directly for read and ref without the whole .c_str() conversion? Could you clarify what you mean by spurious? Do you mean that affine, banded, affine-banded behave inconsistently (and wrong)? Or do you mean that multiple executions of the same binary produces different output? Anyway, I put @rrahn <https://github.com/rrahn> which is the maintainer of the alignment code and should fairly easy find the issue. (Who is on vacation 😄 ) — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#2323 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABTzPuFbb0uckdQAvAgs9Wuyp21Tyzpeks5uHa-FgaJpZM4VG43-> .

[rrahn]

hi, many thanks for reporting this.
I will look into it, but at the moment I am stacked with things that need to be done before we have our workshops in September. So I can't guarantee to fix this by then, but probably when porting the alignment code to SeqAn3 I will fix this as well.