Function find_complex_overlap

DAssembler/DAssembler.cpp:171–301

Find all overlaps with the given focal read and call a consensus sequence

// Find all overlaps with the given focal read and call a consensus sequence
static string find_complex_overlap(const RotatedRead& f,
        const Rotations& r,
        vector<RotatedRead>& rl)
{
    // A vector for tracking all the overlaps, seeded with the initial read
    vector<Overlap> o;
    o.push_back(Overlap(f.seq, 0));

    // The pre-pended string to use to seed the search
    const string& seq1 = '$' + f.seq;

    //Find it in the sorted list - NOTE: if the flank read doesn't correspond
    // to a real read, this iterator will not be used
    Rotations::const_iterator rt = lower_bound(r.begin(), r.end(), seq1);

    /*
    Continue down the sorted list, checking for other matches
     - for real reads (seq1 == rt->seq), continue from the position
    of seq1 in the list
     - otherwise, just start at the beginning
     */
    for(Rotations::const_iterator st = (seq1 == rt->seq) ? rt+1 : r.begin();
            st != r.end(); ++st)
    {
        // Check for an overlap between the two sequences,
        //  allowing for mismatches
        const string& seq2 = st->seq;
        unsigned new_overlap = tier_overlap(seq1, seq2, true);

        // Continue if there's no match
        if (new_overlap == 0
                || new_overlap > opt::max_overlap)
            continue;

        // Add a new overlap object for each appropriate overlap found
        o.push_back(Overlap(
            original_read_from_rotated_read(seq2), new_overlap));
    }

    // Counters for calculating coverage
    // Vector size should be something like "read_length + maximum tier"
    // THIS WILL BREAK WITH LONGER READS
    vector<BaseCount> counts(300);

    // Pretty-print the alignment for verbose only
    if(opt::verbose){
        cerr << endl;
        sort(o.begin(), o.end(), offset_sort);
    }

    // Go through each overlap to count bases offset by the appropriate amount
    for(vector<Overlap>::const_iterator ot = o.begin();
            ot != o.end(); ++ot){

        // Pretty-print each found read aligned with the focal read
        if (opt::verbose)
            cerr << string(ot->offset, ' ') << ot->seq << " t:" << ot->offset;