nw_spliced_aligner32.cpp

ncbi源码

                        V = vAcc;
                        dnr_pos = k0 + jBestDonor[st];
                    }
                }
                else {   // try arbitrary splice
                    TScore vAcc = vBestDonor[st] + m_Wi[st] + m_Wd2;
                    if(vAcc > V) {
                        V = vAcc;
                        dnr_pos = k0 + jBestDonor[st];
                    }
                }
            }
            
            if(dnr_pos != kMax_UI4) {
                type = kTypeIntron | dnr_pos;
            }

            backtrace_matrix[k] = type;

            // detect donor candidates
            if(j < N2 - 2) {
                unsigned char d1 = seq2[j+1], d2 = seq2[j+2];
                Uint1 dnr_mask = 0xF0 & dnr_acc_matrix[(size_t(d1)<<8)|d2];
                for(Uint1 st = 0; st < splice_type_count_32; ++st ) {
                    if( dnr_mask & (0x10 << st) ) {
                        if( d1 == g_nwspl32_donor[st][0] && d2 == g_nwspl32_donor[st][1] ) {
                            if(V > vBestDonor[st]) {
                                jAllDonors[st][jTail[st]] = j;
                                vAllDonors[st][jTail[st]] = V;
                                ++(jTail[st]);
                            }
                        } else {
                            TScore v = V + m_Wd1;
                            if(v > vBestDonor[st]) {
                                jAllDonors[st][jTail[st]] = j;
                                vAllDonors[st][jTail[st]] = v;
                                ++(jTail[st]);
                            }
                        }
                    }
                    else { // both chars distorted
                        TScore v = V + m_Wd2;
                        if(v > vBestDonor[st]) {
                            jAllDonors[st][jTail[st]] = j;
                            vAllDonors[st][jTail[st]] = v;
                            ++(jTail[st]);
                        }
                    }
                }
            }
        }

        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 CSplicedAligner32::x_DoBackTrace ( const Uint4* backtrace_matrix,
                                         size_t N1, size_t N2,
                                         vector<ETranscriptSymbol>* transcript)
{
    transcript->clear();
    transcript->reserve(N1 + N2);

    const Uint4 mask_jump = 0x3FFFFFFF;
    const Uint4 mask_type = ~mask_jump;

    size_t k = N1*N2 - 1;

    while (k != 0) {
        Uint4 Key = backtrace_matrix[k];
        Uint4 type = Key & mask_type;

	if(type == kTypeDiag) {
	    transcript->push_back(eTS_Match);
	    k -= N2 + 1;
	}
	else {
	  Uint4 k2 = (Key & mask_jump);
	  if(k2 >= k) {
	        NCBI_THROW(CAlgoAlignException, eInternal,
                           "Incorrect backtrace jump detected");
	  }

	  ETranscriptSymbol ts;
	  Uint4 decr;
	  if(type == kTypeIntron) {
	    ts = eTS_Intron;
	    decr = 1;
	  }
	  else {
	    Uint4 kdel = k / N2, k2del = k2 / N2;
	    Uint4 kmod = k % N2, k2mod = k2 % N2;
	    if(kdel == k2del) {
	      ts = eTS_Insert;
	      decr = 1;
	    }
	    else if(kmod == k2mod) {
	      ts = eTS_Delete;
	      decr = N2;
	    }
	    else {
	      // ts = eTS_DiagSpace;
	      NCBI_THROW(CAlgoAlignException, eInternal,
                         "Diag spaces not yet supported");
	    }
	  }

	  for(; k > k2; k -= decr) {
	    transcript->push_back(ts);
	  }
	  if(k != k2) {
	      NCBI_THROW(CAlgoAlignException, eInternal,
                         "Backtrace jump failed");
	  }
	}
    } // while
}



CNWAligner::TScore CSplicedAligner32::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 k = 0; k < splice_type_count_32; ++k) {
                    if(*p2 == g_nwspl32_donor[k][0] && *(p2 + 1) == g_nwspl32_donor[k][1]) {
                        score += m_Wi[k];
                        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_aligner32.cpp,v $
 * Revision 1000.3  2004/06/01 18:05:06  gouriano
 * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.15
 *
 * Revision 1.15  2004/05/21 21:41:02  gorelenk
 * Added PCH ncbi_pch.hpp
 *
 * Revision 1.14  2004/05/18 21:43:40  kapustin
 * Code cleanup
 *
 * Revision 1.13  2004/05/17 14:50:57  kapustin
 * Add/remove/rearrange some includes and object declarations
 *
 * Revision 1.12  2004/04/23 14:39:47  kapustin
 * Add Splign library and other changes
 *
 * Revision 1.10  2003/11/20 17:54:23  kapustin
 * Alternative conventional splice penalty adjusted
 *
 * Revision 1.9  2003/10/31 19:40:13  kapustin
 * Get rid of some WS and GCC complains
 *
 * Revision 1.8  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.7  2003/10/14 19:31:52  kapustin
 * Use one flag for all gap types
 *
 * Revision 1.6  2003/09/30 19:50:04  kapustin
 * Make use of standard score matrix interface
 *
 * Revision 1.5  2003/09/26 14:43:18  kapustin
 * Remove exception specifications
 *
 * Revision 1.4  2003/09/10 20:25:21  kapustin
 * Use unsigned char for score matrix index
 *
 * Revision 1.3  2003/09/10 19:13:10  kapustin
 * Change backtrace type constants to allow larger jump values
 *
 * 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
 *
 * ===========================================================================
 */