mm_aligner.cpp

ncbi源码

        unsigned char ci = seq1[i];

        TScore wg2 = m_Wg, ws2 = m_Ws;

        for (j = 1; j < N2; ++j) {

            G = pV[j] + sm[ci][(unsigned char)seq2[j]];
            pV[j] = V;

            n0 = V + m_Wg;
            if(E >= n0)
                E += m_Ws;      // continue the gap
            else
                E = n0 + m_Ws;  // open a new gap

            if(j == N2 - 1 && bFreeGapRight2)
                wg2 = ws2 = 0;

            n0 = rowV[j] + wg2;
            if (rowF[j] > n0)
                rowF[j] += ws2;
            else
                rowF[j] = n0 + ws2;

            V = (E >= rowF[j])? (E >= G? E: G): (rowF[j] >= G? rowF[j]: G);
        }
        pV[j] = V;

        if( m_prg_callback && i % prg_rep_rate == 0 ) {
#ifdef NCBI_THREADS
            CFastMutexGuard guard (progress_mutex);
#endif
            m_prg_info.m_iter_done += prg_rep_increment;
            if( m_terminate = m_prg_callback(&m_prg_info) ) {
                break;
            }
        }
    }

    // the last row (i == N1 - 1)
    if(!m_terminate) {

        vG[0] = vE[0] = E = kInfMinus;
        vF[0] = V = V0 += ws;
        trace[0] = kMaskFc;
        unsigned char ci = seq1[i];

        TScore wg2 = m_Wg, ws2 = m_Ws;

        unsigned char tracer;
        for (j = 1; j < N2; ++j) {

            vG[j] = G = pV[j] + sm[ci][(unsigned char)seq2[j]];
            pV[j] = V;

            n0 = V + m_Wg;
            if(E >= n0) {
                E += m_Ws;      // continue the gap
                tracer = kMaskEc;
            }
            else {
                E = n0 + m_Ws;  // open a new gap
                tracer = 0;
            }
            vE[j] = E;

            if(j == N2 - 1 && bFreeGapRight2) {
                wg2 = ws2 = 0;
            }

            n0 = rowV[j] + wg2;
            if(rowF[j] >= n0) {
                rowF[j] += ws2;
                tracer |= kMaskFc;
            }
            else {
                rowF[j] = n0 + ws2;
            }
            vF[j] = rowF[j];

            if (E >= rowF[j]) {
                if(E >= G) {
                    V = E;
                    tracer |= kMaskE;
                }
                else {
                    V = G;
                    tracer |= kMaskD;
                }
            } else {
                if(rowF[j] >= G) {
                    V = rowF[j];
                }
                else {
                    V = G;
                    tracer |= kMaskD;
                }
            }
            trace[j] = tracer;
        }
    }

    if( m_prg_callback ) {
#ifdef NCBI_THREADS
        CFastMutexGuard guard (progress_mutex);
#endif
        m_prg_info.m_iter_done += (i % prg_rep_rate)*N2;
        m_terminate = m_prg_callback(&m_prg_info);
    }
}

 
void CMMAligner::x_RunBtm(const SCoordRect& rect,
             vector<TScore>& vE, vector<TScore>& vF, vector<TScore>& vG,
             vector<unsigned char>& trace, bool rb) const
{
    if( m_terminate ) {
        return;
    }

    const size_t dim1 = rect.i2 - rect.i1 + 1;
    const size_t dim2 = rect.j2 - rect.j1 + 1;
    const size_t N1   = dim1 + 1;
    const size_t N2   = dim2 + 1;

    vector<TScore> stl_rowV (N2), stl_rowF (N2);
    TScore* rowV    = &stl_rowV [0];
    TScore* rowF    = &stl_rowF [0];

    TScore* pV = rowV + 1;

    const char* seq1 = m_Seq1 + rect.i1;
    const char* seq2 = m_Seq2 + rect.j1;

    const TNCBIScore (*sm) [NCBI_FSM_DIM] = m_ScoreMatrix.s;

    bool bFreeGapRight1  = m_esf_R1 && rect.i2 == m_SeqLen1 - 1;
    bool bFreeGapRight2  = m_esf_R2 && rect.j2 == m_SeqLen2 - 1;
    bool bFreeGapLeft2  =  m_esf_L2 && rect.j1 == 0;

    // progress reporting

    const size_t prg_rep_rate = 100;
    const size_t prg_rep_increment = prg_rep_rate*N2;

    // bottom row

    TScore wg = bFreeGapRight1? 0: m_Wg;
    TScore ws = bFreeGapRight1? 0: m_Ws;

    rowV[N2 - 1] = wg;
    int i, j;
    for (j = N2 - 2; j >= 0; --j) {
        rowV[j] = pV[j] + ws;
        rowF[j] = kInfMinus;
    }
    rowV[N2 - 1] = 0;

    // recurrences

    wg = bFreeGapRight2? 0: m_Wg;
    ws = bFreeGapRight2? 0: m_Ws;

    TScore V  = 0;
    TScore V0 = rb? 0: wg;
    TScore E, G, n0;

    for(i = N1 - 2;  i > 0;  --i) {
        
        V = V0 += ws;
        E = kInfMinus;
        unsigned char ci = seq1[i];

        TScore wg2 = m_Wg, ws2 = m_Ws;

        for (j = N2 - 2; j >= 0; --j) {

            G = pV[j] + sm[ci][(unsigned char)seq2[j]];
            pV[j] = V;

            n0 = V + m_Wg;
            if(E >= n0)
                E += m_Ws;      // continue the gap
            else
                E = n0 + m_Ws;  // open a new gap

            if(j == 0 && bFreeGapLeft2) {
                wg2 = ws2 = 0;
            }

            n0 = rowV[j] + wg2;
            if (rowF[j] > n0)
                rowF[j] += ws2;
            else
                rowF[j] = n0 + ws2;

            V = (E >= rowF[j])? (E >= G? E: G): (rowF[j] >= G? rowF[j]: G);
        }
        pV[j] = V;

        if( m_prg_callback && (N1 - i) % prg_rep_rate == 0 ) {
#ifdef NCBI_THREADS
            CFastMutexGuard guard (progress_mutex);
#endif
            m_prg_info.m_iter_done += prg_rep_increment;
            if( m_terminate = m_prg_callback(&m_prg_info) ) {
                break;
            }
        }
    }

    // the top row (i == 0)
    if(!m_terminate) {

        vF[N2-1] = V = V0 += ws;
        vG[N2-1] = vE[N2-1] = E = kInfMinus;
        trace[N2-1] = kMaskFc;
        unsigned char ci = seq1[i];

        TScore wg2 = m_Wg, ws2 = m_Ws;

        unsigned char tracer;
        for (j = N2 - 2; j >= 0; --j) {

            vG[j] = G = pV[j] + sm[ci][(unsigned char)seq2[j]];
            pV[j] = V;

            n0 = V + m_Wg;
            if(E >= n0) {
                E += m_Ws;      // continue the gap
                tracer = kMaskEc;
            }
            else {
                E = n0 + m_Ws;  // open a new gap
                tracer = 0;
            }
            vE[j] = E;

            if(j == 0 && bFreeGapLeft2) {
                wg2 = ws2 = 0;
            }

            n0 = rowV[j] + wg2;
            if(rowF[j] >= n0) {
                rowF[j] += ws2;
                tracer |= kMaskFc;
            }
            else {
                rowF[j] = n0 + ws2;
            }
            vF[j] = rowF[j];

            if (E >= rowF[j]) {
                if(E >= G) {
                    V = E;
                    tracer |= kMaskE;
                }
                else {
                    V = G;
                    tracer |= kMaskD;
                }
            } else {
                if(rowF[j] >= G) {
                    V = rowF[j];
                }
                else {
                    V = G;
                    tracer |= kMaskD;
                }
            }
            trace[j] = tracer;
        }
    }

    if( m_prg_callback ) {
#ifdef NCBI_THREADS
        CFastMutexGuard guard (progress_mutex);
#endif
        m_prg_info.m_iter_done += (N1 - i) % prg_rep_rate;
        m_terminate = m_prg_callback(&m_prg_info);
    }
}


CNWAligner::TScore CMMAligner::x_RunTerm(const SCoordRect& rect,
                                         bool left_top, bool right_bottom,
                                         list<ETranscriptSymbol>& subpath)
{
    if( m_terminate ) {
        return 0;
    }

    const size_t N1 = rect.i2 - rect.i1 + 2;
    const size_t N2 = rect.j2 - rect.j1 + 2;

    vector<TScore> stl_rowV (N2), stl_rowF (N2);
    TScore* rowV    = &stl_rowV [0];
    TScore* rowF    = &stl_rowF [0];

    // index calculation: [i,j] = i*n2 + j
    vector<unsigned char> stl_bm (N1*N2);
    unsigned char* backtrace = &stl_bm[0];

    TScore* pV = rowV - 1;

    const char* seq1 = m_Seq1 + rect.i1 - 1;
    const char* seq2 = m_Seq2 + rect.j1 - 1;

    const TNCBIScore (*sm) [NCBI_FSM_DIM] = m_ScoreMatrix.s;

    bool bFreeGapLeft1  = m_esf_L1 && rect.i1 == 0;
    bool bFreeGapRight1 = m_esf_R1 && rect.i2 == m_SeqLen1 - 1;
    bool bFreeGapLeft2  = m_esf_L2 && rect.j1 == 0;
    bool bFreeGapRight2 = m_esf_R2 && rect.j2 == m_SeqLen2 - 1;

    TScore wgleft1   = bFreeGapLeft1? 0: m_Wg;
    TScore wsleft1   = bFreeGapLeft1? 0: m_Ws;
    TScore wg1 = m_Wg, ws1 = m_Ws;

    // first row
    size_t k;
    {
        rowV[0] = wgleft1;
        for (k = 1; k < N2; k++) {
            rowV[k] = pV[k] + wsleft1;
            rowF[k] = kInfMinus;
            backtrace[k] = kMaskE | kMaskEc;
        }
        rowV[0] = 0;
    }

    // recurrences
    TScore wgleft2   = bFreeGapLeft2? 0: m_Wg;
    TScore wsleft2   = bFreeGapLeft2? 0: m_Ws;
    TScore V  = 0;
    TScore V0 = left_top? 0: wgleft2;
    TScore E, G, n0;
    unsigned char tracer;

    size_t i, j;
    for(i = 1;  i < N1;  ++i) {
        
        V = V0 += wsleft2;
        E = kInfMinus;
        backtrace[k++] = kMaskFc;
        unsigned char ci = seq1[i];

        if(i == N1 - 1 && bFreeGapRight1) {
                wg1 = ws1 = 0;
        }

        TScore wg2 = m_Wg, ws2 = m_Ws;

        for (j = 1; j < N2; ++j, ++k) {

            G = pV[j] + sm[ci][(unsigned char)seq2[j]];
            pV[j] = V;

            n0 = V + wg1;
            if(E >= n0) {
                E += ws1;      // continue the gap
                tracer = kMaskEc;
            }
            else {
                E = n0 + ws1;  // open a new gap
                tracer = 0;
            }

            if(j == N2 - 1 && bFreeGapRight2) {
                wg2 = ws2 = 0;
            }
            n0 = rowV[j] + ((right_bottom && j == N2 - 1)? 0: wg2);
            if(rowF[j] >= n0) {
                rowF[j] += ws2;
                tracer |= kMaskFc;
            }
            else {
                rowF[j] = n0 + ws2;
            }

            if (E >= rowF[j]) {
                if(E >= G) {
                    V = E;
                    tracer |= kMaskE;
                }
                else {
                    V = G;
                    tracer |= kMaskD;
                }
            } else {
                if(rowF[j] >= G) {
                    V = rowF[j];
                }
                else {
                    V = G;
                    tracer |= kMaskD;
                }
            }
            backtrace[k] = tracer;
        }

        pV[j] = V;
    }

    // fill the subpath
    subpath.clear();
    
    // run backtrace
    k = N1*N2 - 1;
    while (k != 0) {
        unsigned char Key = backtrace[k];
        if (Key & kMaskD) {
            subpath.push_front(eTS_Match);
            k -= N2 + 1;
        }
        else if (Key & kMaskE) {
            subpath.push_front(eTS_Insert); --k;
            while(k > 0 && (Key & kMaskEc)) {
                subpath.push_front(eTS_Insert);
                Key = backtrace[k--];
            }
        }
        else {
            subpath.push_front(eTS_Delete);
            k -= N2;
            while(k > 0 && (Key & kMaskFc)) {
                subpath.push_front(eTS_Delete);
                Key = backtrace[k];
                k -= N2;
            }
        }
    }

    return V;
}


bool CMMAligner::x_CheckMemoryLimit()
{
    return true;
}



END_NCBI_SCOPE

/*
 * ===========================================================================
 * $Log: mm_aligner.cpp,v $
 * Revision 1000.1  2004/06/01 18:04:46  gouriano
 * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.20
 *
 * Revision 1.20  2004/05/21 21:41:01  gorelenk
 * Added PCH ncbi_pch.hpp
 *
 * Revision 1.19  2004/05/18 21:43:40  kapustin
 * Code cleanup
 *
 * Revision 1.18  2004/05/17 14:50:56  kapustin
 * Add/remove/rearrange some includes and object declarations
 *
 * Revision 1.17  2003/09/30 19:50:04  kapustin
 * Make use of standard score matrix interface
 *
 * Revision 1.16  2003/09/26 14:43:18  kapustin
 * Remove exception specifications
 *
 * Revision 1.15  2003/09/04 16:07:38  kapustin
 * Use STL vectors for exception-safe dynamic arrays and matrices
 *
 * Revision 1.14  2003/09/02 22:36:57  kapustin
 * Adjust for transcript symbol name changes
 *
 * Revision 1.13  2003/06/17 17:20:44  kapustin
 * CNWAlignerException -> CAlgoAlignException
 *
 * Revision 1.12  2003/06/17 14:51:04  dicuccio
 * Fixed after algo/ rearragnement
 *
 * Revision 1.11  2003/06/03 18:19:06  rsmith
 * Move static mutex initialization to file from function scope,
 * because of MW compiler choking (wrongly) over complex initialization.
 *
 * Revision 1.10  2003/06/02 14:04:49  kapustin
 * Progress indication-related updates
 *
 * Revision 1.9  2003/05/23 18:26:38  kapustin
 * Use weak comparisons in core recurrences.
 *
 * Revision 1.8  2003/04/14 19:01:49  kapustin
 * Run() -> x_Run()
 *
 * Revision 1.7  2003/03/18 15:15:51  kapustin
 * Implement virtual memory checking function. Allow separate free
 * end gap specification
 *
 * Revision 1.6  2003/03/17 15:30:57  kapustin
 * Support end-space free alignments
 *
 * Revision 1.5  2003/01/30 20:32:06  kapustin
 * Call x_LoadScoringMatrix() from the base class constructor.
 *
 * Revision 1.4  2003/01/28 12:37:40  kapustin
 * Move m_score to the base class
 *
 * Revision 1.3  2003/01/22 13:40:09  kapustin
 * Implement multithreaded algorithm
 *
 * ===========================================================================
 */