nw_spliced_aligner16.cpp

ncbi源码

            }

            if(tracer_dnr != 0xFFFF) {
                backtrace_matrix[dnr_pos] |= tracer_dnr;
                tracer |= tracer_acc;
            }

            backtrace_matrix[k] = tracer;

            // detect donor candidate
            if(j < N2 - 2) {
                for(unsigned char st = 0; st < g_topidx; ++st) {
                    
                    if( seq2[j+1] == g_nwspl_donor[st][0] &&
                        seq2[j+2] == g_nwspl_donor[st][1] &&
                        V > vBestDonor[st]  ) {
                        
                        jAllDonors[st][jTail[st]] = j;
                        vAllDonors[st][jTail[st]] = V;
                        ++(jTail[st]);
                    }
                }
            }
            
            // detect new best value
            if(V > V_max) {
                jAllDonors[g_topidx][jTail[g_topidx]] = j;
                vAllDonors[g_topidx][jTail[g_topidx]] = V_max = V;
                ++(jTail[3]);
            }
        }

        pV[j] = V;

        if(i == 0) {
            V0 = wgleft2;
            wg1 = m_Wg;
            ws1 = m_Ws;
        }

    }

    try {
        x_DoBackTrace(backtrace_matrix, N1, N2, transcript);
    }
    catch(exception&) { // GCC hack
        throw;
    }
    
    return V;
}


// perform backtrace step;
void CSplicedAligner16::x_DoBackTrace ( const Uint2* backtrace_matrix,
                                         size_t N1, size_t N2,
                                         vector<ETranscriptSymbol>* transcript)
{
    transcript->clear();
    transcript->reserve(N1 + N2);

    size_t k = N1*N2 - 1;
    while (k != 0) {
        Uint2 Key = backtrace_matrix[k];
        while(Key & 0x0070) {  // acceptor

            unsigned char acc_type = (Key & 0x0070) >> 4;
            Uint2 dnr_mask = 0x0040 << acc_type;
            ETranscriptSymbol ets = eTS_Intron;
            do {
                transcript->push_back(ets);
                Key = backtrace_matrix[--k];
            }
            while(!(Key & dnr_mask));

            if(Key & 0x0070) {
                NCBI_THROW(CAlgoAlignException,
                           eInternal,
                           "Adjacent splices encountered during backtrace");
            }
        }

        if(k == 0) continue;
        
        if (Key & kMaskD) {
            transcript->push_back(eTS_Match);  // or eTS_Replace
            k -= N2 + 1;
        }
        else if (Key & kMaskE) {
            transcript->push_back(eTS_Insert); --k;
            while(k > 0 && (Key & kMaskEc)) {
                transcript->push_back(eTS_Insert);
                Key = backtrace_matrix[k--];
            }
        }
        else {
            transcript->push_back(eTS_Delete);
            k -= N2;
            while(k > 0 && (Key & kMaskFc)) {
                transcript->push_back(eTS_Delete);
                Key = backtrace_matrix[k];
                k -= N2;
            }
        }
    }
}


CNWAligner::TScore CSplicedAligner16::x_ScoreByTranscript() const
{
    const size_t dim = m_Transcript.size();
    if(dim == 0) return 0;

    TScore score = 0;

    const char* p1 = m_Seq1;
    const char* p2 = m_Seq2;

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

    char state1;   // 0 = normal, 1 = gap, 2 = intron
    char state2;   // 0 = normal, 1 = gap

    switch( m_Transcript[dim - 1] ) {

    case eTS_Match:
        state1 = state2 = 0;
        break;

    case eTS_Insert:
        state1 = 1; state2 = 0; score += m_Wg;
        break;

    case eTS_Intron:
        state1 = 0; state2 = 0;
        break; // intron flag set later

    case eTS_Delete:
        state1 = 0; state2 = 1; score += m_Wg;
        break;

    default: {
        NCBI_THROW(CAlgoAlignException, eInternal,
                   "Invalid transcript symbol");
        }
    }

    for(int i = dim - 1; i >= 0; --i) {

        unsigned char c1 = m_Seq1? *p1: 0;
        unsigned char c2 = m_Seq2? *p2: 0;
        switch(m_Transcript[i]) {

        case eTS_Match: {
            state1 = state2 = 0;
            score += sm[c1][c2];
            ++p1; ++p2;
        }
        break;

        case eTS_Insert: {
            if(state1 != 1) score += m_Wg;
            state1 = 1; state2 = 0;
            score += m_Ws;
            ++p2;
        }
        break;

        case eTS_Delete: {
            if(state2 != 1) score += m_Wg;
            state1 = 0; state2 = 1;
            score += m_Ws;
            ++p1;
        }
        break;

        case eTS_Intron: {

            if(state1 != 2) {
                for(unsigned char i = 0; i < splice_type_count_16; ++i) {
                    if(*p2 == g_nwspl_donor[i][0] &&
                       *(p2 + 1) == g_nwspl_donor[i][1] || i == g_topidx) {

                        score += m_Wi[i];
                        break;
                    }
                }
            }
            state1 = 2; state2 = 0;
            ++p2;
        }
        break;

        default: {
        NCBI_THROW(CAlgoAlignException, eInternal,
                   "Invalid transcript symbol");
        }
        }
    }

    if(m_esf_L1) {
        size_t g = 0;
        for(int i = dim - 1; i >= 0; --i) {
            if(m_Transcript[i] == eTS_Insert) ++g; else break;
        }
        if(g > 0) {
            score -= (m_Wg + g*m_Ws);
        }
    }

    if(m_esf_L2) {
        size_t g = 0;
        for(int i = dim - 1; i >= 0; --i) {
            if(m_Transcript[i] == eTS_Delete) ++g; else break;
        }
        if(g > 0) {
            score -= (m_Wg + g*m_Ws);
        }
    }

    if(m_esf_R1) {
        size_t g = 0;
        for(size_t i = 0; i < dim; ++i) {
            if(m_Transcript[i] == eTS_Insert) ++g; else break;
        }
        if(g > 0) {
            score -= (m_Wg + g*m_Ws);
        }
    }

    if(m_esf_R2) {
        size_t g = 0;
        for(size_t i = 0; i < dim; ++i) {
            if(m_Transcript[i] == eTS_Delete) ++g; else break;
        }
        if(g > 0) {
            score -= (m_Wg + g*m_Ws);
        }
    }

    return score;
}


END_NCBI_SCOPE



/*
 * ===========================================================================
 * $Log: nw_spliced_aligner16.cpp,v $
 * Revision 1000.3  2004/06/01 18:05:00  gouriano
 * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.13
 *
 * Revision 1.13  2004/05/21 21:41:02  gorelenk
 * Added PCH ncbi_pch.hpp
 *
 * Revision 1.12  2004/05/18 21:43:40  kapustin
 * Code cleanup
 *
 * Revision 1.11  2004/05/17 14:50:56  kapustin
 * Add/remove/rearrange some includes and object declarations
 *
 * Revision 1.10  2004/04/23 14:39:47  kapustin
 * Add Splign library and other changes
 *
 * Revision 1.8  2003/11/20 17:54:23  kapustin
 * Alternative conventional splice penalty adjusted
 *
 * Revision 1.7  2003/10/31 19:40:13  kapustin
 * Get rid of some WS and GCC complains
 *
 * Revision 1.6  2003/10/27 21:00:17  kapustin
 * Set intron penalty defaults differently for 16- and 32-bit versions
 * according to the expected quality of sequences those variants are
 * supposed to be used with.
 *
 * Revision 1.5  2003/10/14 19:29:24  kapustin
 * Dismiss static keyword as a local-to-compilation-unit flag. Use longer
 * name since unnamed namespaces are not everywhere supported
 *
 * Revision 1.4  2003/09/30 19:50:04  kapustin
 * Make use of standard score matrix interface
 *
 * Revision 1.3  2003/09/26 14:43:18  kapustin
 * Remove exception specifications
 *
 * Revision 1.2  2003/09/04 16:07:38  kapustin
 * Use STL vectors for exception-safe dynamic arrays and matrices
 *
 * Revision 1.1  2003/09/02 22:34:49  kapustin
 * Initial revision
 *
 * ===========================================================================
 */