[Seqan-dev] Fatal Crash / Failed Assert in BandedAlignment

[Brett Bowman <bnbowman@gmail.com>]

I'm using Banded aligner to find potentially low-quality but important spike-in control alignments from single-molecule sequencing data (PacBio, Oxford Nano).  When these potential alignments turn out to be junk, bandedAlignment is sometimes either failing on an assert (Debug) or just out-right seg-faulting (Release), instead of getting the expected junk alignment or graceful failure I would expect. 

As it stands, I've implemented some filters to toss these low-quality seed-chains out before alignment, but it would be preferable if I didn't have to worry about the aligner crashing at all...

Sample Error Message:

"""

/home/UNIXHOME/bbowman/p4/depot/software/smrtanalysis/bioinformatics/third-party/cpp/seqan-library-2.0.0/include/seqan/seeds/banded_chain_alignment_traceback.h:195 Assertion failed : isGlued == true was: 0 != 1 (Fatal error while trying to connect trace backs: No glue point available!)

"""

Sample Test Data/Code:

"""

#include <seqan/seeds.h>

#include <seqan/journaled_set/score_biaffine.h>

// Type-defs used

typedef seqan::Dna TAlphabet;

typedef seqan::String<TAlphabet> TString;

typedef seqan::Seed<seqan::Simple> TSeed;

typedef seqan::String<TSeed> TSeedChain;

typedef seqan::AlignConfig<true, true, true, true> TAlignConfig;

typedef seqan::Align<TString, seqan::ArrayGaps> TAlign;

typedef seqan::Score<long, seqan::BiAffine> TScoringScheme;

// Variables used

const TString query = "AAAGAGAGAGAGAGAGAGGGAAAAGAAGGAAAGGGAAAAGGAGAAGGGAGGGAAAAAGAAAGAAAGAAGAAAAGAGAGGGAAAAGAGGAAGAAGAAAGGGAGGGGAAGAGGAAGAGGGAGAGGAGAGAAGGAGAGGAGGAAGAAAAAAAAAAGAGAAAAGGGGAGGAGAGGAGAAAGGAAGAAGGGGAAAAGAGGGGAGAGGAAGGAGAAGAGGAAAGAAAAAAAAGGAGAGAAAGGGAAGGAGAAAAGAGGAAGGAGGAGAGAAGAAAGAAAGAGGAGAAAAGAGAAGGAAAAAGAAAGGAAGGAGGGAGAGGAGAGAGGAGAAAAGAAAAGGGAAAAGAAAAAGGAAAGAAAAAGAAAGGAAAAAGAAAGAGAAGAAAAGGAAAAAGAGAAGGGAGGAAAAGAGGAAGAGGAGAAAAAAAAAGGGAGAAAAAAGAGAAAAGGGAGGAAGAAAAGAAAAAAGGAGGGAGGAGAAAGAAAGAGAAAGAAAGAGGGGAGAAGAAGAGGAAGGAAGGGAGAGAAAGAGGAAAAAAAGGGAAAGAGGAGAGGAAAGAGGAGAGAAAGAGGAAAGAAGAAAAGAGAAAAAAAAAAAAGAAGAAGAGAAGAAAAAAAGAAAAAGGGAAAAAGAAGAGGGAAAGGAAAAGAGAAGAGAAAGGGAGAGGAAAAGAGAGAAAGAGAAAAGAGGAAAAAAAAGAGAGGAAGGAGGAAAAGAAAAAAAGGAA";

const TString ref   = "AAAGAGAGAGAGAAGGCTGGGCAGGCTATGCACCCTGGTCCAGGTCAAAAGCTGCGGAACCCGCTAGCGGCCATCTTGGCCACTAGGGGTCCCGCAGATTCATATTGTCGTCTAGCATGCACAATGCTGCAAACCCAGCTTGCAATGCCCACAGCAAGCGGCCAATCTTTACGCCACGTTGAATTGTTTATTACCTGTGACTGGCTATGGCTTGCAACGCCACTCGTAAAACTAGTACTTTGCGGTTAGGGGAAGTAGACAAACCCATTACTCCACTTCCCGGAAGTTCAACTCATTCCAACACGAAATAAAAGTAAACTCAACACCCCAAGCAGGCTATGTGGGGGGGTGATAGGGGTGGATTCTATTTCCTATCCCATCCCCTAGGATCTCAATTAAGTTACTAGCGAGTTAAATGTCTGTAGCGATCCCGTCAGTCCTATCGCGCGCATCAAGACCTGGTTGGTTGAGCGTGCAGTAGATCATCGATAAGCTGCGAGTTAGGTCATCCCAGACCGCATCTGGCGCCTAAACGTTCAGTGGTAGCTAAGGCGTCACCTTCGACTGTCTAAAGGCAATATGTCGTCCTTAGCTCCAAGTCCCTAGCAAGCGTGTCGGGTCTCTCTCAACAACAACAACGGAGGAGGAGGAAAAGAGAGAGACCCGACACGCTTGCTAGGGACTTGGAGCTAAGGACGACATATTGCCTT";

const TScoringScheme scoringScheme = TScoringScheme(5, -6, -4, -4, -5, -5);

const TAlignConfig alignConfig     = TAlignConfig();

//Seed String of Type: Shift 1

//Seed #1 Seed<Simple, TConfig>(1, 1, 13, 13, lower diag = 0, upper diag = 0)

//Seed #2 Seed<Simple, TConfig>(660, 648, 672, 660, lower diag = 12, upper diag = 12)

int main(int argc, char** argv)

{

    // Seed chain

    TSeedChain chain;

    seqan::appendValue(chain, TSeed(1,     1, 13));

    seqan::appendValue(chain, TSeed(660, 648, 12));

    // Alignment

    TAlign align;

    seqan::resize(seqan::rows(align), 2);

    seqan::assignSource(seqan::row(align, 0), query);

    seqan::assignSource(seqan::row(align, 1),   ref);

    bandedChainAlignment(align, chain, scoringScheme, alignConfig, 15);

    return 0;

}

"""

Sincerely,

Brett