开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 生物技术 > 查看源码
nw_aligner.cpp
资源名称：ncbi_cxx_win--Jul_08_2004.zip　[点击查看]
上传用户：yhdzpy8989
上传日期：2007-06-13
资源大小：13604k
文件大小：31k
源码类别：

生物技术

开发平台：

C/C++

nw_aligner.cpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: nw_aligner.cpp,v $
  4.  * PRODUCTION Revision 1000.3  2004/06/01 18:04:52  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.48
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /* $Id: nw_aligner.cpp,v 1000.3 2004/06/01 18:04:52 gouriano Exp $
  11.  * ===========================================================================
  12.  *
  13.  *                            PUBLIC DOMAIN NOTICE                          
  14.  *               National Center for Biotechnology Information
  15.  *                                                                          
  16.  *  This software/database is a "United States Government Work" under the   
  17.  *  terms of the United States Copyright Act.  It was written as part of    
  18.  *  the author's official duties as a United States Government employee and 
  19.  *  thus cannot be copyrighted.  This software/database is freely available 
  20.  *  to the public for use. The National Library of Medicine and the U.S.    
  21.  *  Government have not placed any restriction on its use or reproduction.  
  22.  *                                                                          
  23.  *  Although all reasonable efforts have been taken to ensure the accuracy  
  24.  *  and reliability of the software and data, the NLM and the U.S.           
  25.  *  Government do not and cannot warrant the performance or results that    
  26.  *  may be obtained by using this software or data. The NLM and the U.S.    
  27.  *  Government disclaim all warranties, express or implied, including       
  28.  *  warranties of performance, merchantability or fitness for any particular
  29.  *  purpose.                                                                
  30.  *                                                                          
  31.  *  Please cite the author in any work or product based on this material.   
  32.  *
  33.  * ===========================================================================
  34.  *
  35.  * Authors:  Yuri Kapustin, Boris Kiryutin
  36.  *
  37.  * File Description:  CNWAligner implementation
  38.  *                   
  39.  * ===========================================================================
  40.  *
  41.  */
  44. #include <ncbi_pch.hpp>
  45. #include <algo/align/nw_aligner.hpp>
  46. #include <algo/align/align_exception.hpp>
  49. BEGIN_NCBI_SCOPE
  51. // IUPACna alphabet (default)
  52. const char g_nwaligner_nucleotides [] = "AGTCBDHKMNRSVWY";
  55. CNWAligner::CNWAligner()
  56.     : m_Wm(GetDefaultWm()),
  57.       m_Wms(GetDefaultWms()),
  58.       m_Wg(GetDefaultWg()),
  59.       m_Ws(GetDefaultWs()),
  60.       m_esf_L1(false), m_esf_R1(false), m_esf_L2(false), m_esf_R2(false),
  61.       m_abc(g_nwaligner_nucleotides),
  62.       m_prg_callback(0),
  63.       m_Seq1(0), m_SeqLen1(0),
  64.       m_Seq2(0), m_SeqLen2(0),
  65.       m_score(kInfMinus)
  66. {
  67.     SetScoreMatrix(0);
  68. }
  71. CNWAligner::CNWAligner( const char* seq1, size_t len1,
  72.                         const char* seq2, size_t len2,
  73.                         const SNCBIPackedScoreMatrix* scoremat )
  75.     : m_Wm(GetDefaultWm()),
  76.       m_Wms(GetDefaultWms()),
  77.       m_Wg(GetDefaultWg()),
  78.       m_Ws(GetDefaultWs()),
  79.       m_esf_L1(false), m_esf_R1(false), m_esf_L2(false), m_esf_R2(false),
  80.       m_abc(g_nwaligner_nucleotides),
  81.       m_prg_callback(0),
  82.       m_Seq1(seq1), m_SeqLen1(len1),
  83.       m_Seq2(seq2), m_SeqLen2(len2),
  84.       m_score(kInfMinus)
  85. {
  86.     SetScoreMatrix(scoremat);
  87.     SetSequences(seq1, len1, seq2, len2);
  88. }
  91. void CNWAligner::SetSequences(const char* seq1, size_t len1,
  92.       const char* seq2, size_t len2,
  93.       bool verify)
  94. {
  95.     if(!seq1 || !seq2) {
  96.         NCBI_THROW(CAlgoAlignException, eBadParameter,
  97.                    "NULL sequence pointer(s) passed");
  98.     }
  100.     if(verify) {
  101.         size_t iErrPos1 = x_CheckSequence(seq1, len1);
  102. if(iErrPos1 < len1) {
  103.     ostrstream oss;
  104.     oss << "The first sequence is inconsistent with the current "
  105. << "scoring matrix type. Symbol " << seq1[iErrPos1] << " at "
  106. << iErrPos1;
  107.     string message = CNcbiOstrstreamToString(oss);
  108.     NCBI_THROW(CAlgoAlignException, eInvalidCharacter, message);
  109. }
  111. size_t iErrPos2 = x_CheckSequence(seq2, len2);
  112. if(iErrPos2 < len2) {
  113.     ostrstream oss;
  114.     oss << "The second sequence is inconsistent with the current "
  115. << "scoring matrix type. Symbol " << seq2[iErrPos2] << " at "
  116. << iErrPos2;
  117.     string message = CNcbiOstrstreamToString(oss);
  118.     NCBI_THROW(CAlgoAlignException, eInvalidCharacter, message);
  119. }
  120.     }
  121.     m_Seq1 = seq1;
  122.     m_SeqLen1 = len1;
  123.     m_Seq2 = seq2;
  124.     m_SeqLen2 = len2;
  125.     m_Transcript.clear();
  126. }
  130. void CNWAligner::SetEndSpaceFree(bool Left1, bool Right1,
  131.                                  bool Left2, bool Right2)
  132. {
  133.     m_esf_L1 = Left1;
  134.     m_esf_R1 = Right1;
  135.     m_esf_L2 = Left2;
  136.     m_esf_R2 = Right2;
  137. }
  139. // evaluate score for each possible alignment;
  140. // fill out backtrace matrix
  141. // bit coding (four bits per value): D E Ec Fc
  142. // D:  1 if diagonal; 0 - otherwise
  143. // E:  1 if space in 1st sequence; 0 if space in 2nd sequence
  144. // Ec: 1 if gap in 1st sequence was extended; 0 if it is was opened
  145. // Fc: 1 if gap in 2nd sequence was extended; 0 if it is was opened
  146. //
  148. const unsigned char kMaskFc  = 0x0001;
  149. const unsigned char kMaskEc  = 0x0002;
  150. const unsigned char kMaskE   = 0x0004;
  151. const unsigned char kMaskD   = 0x0008;
  153. CNWAligner::TScore CNWAligner::x_Align(const char* seg1, size_t len1,
  154.                                        const char* seg2, size_t len2,
  155.                                        vector<ETranscriptSymbol>* transcript)
  156. {
  157.     const size_t N1 = len1 + 1;
  158.     const size_t N2 = len2 + 1;
  160.     vector<TScore> stl_rowV (N2), stl_rowF(N2);
  162.     TScore* rowV    = &stl_rowV[0];
  163.     TScore* rowF    = &stl_rowF[0];
  165.     TScore* pV = rowV - 1;
  167.     const char* seq1 = seg1 - 1;
  168.     const char* seq2 = seg2 - 1;
  170.     const TNCBIScore (*sm) [NCBI_FSM_DIM] = m_ScoreMatrix.s;
  172.     m_terminate = false;
  174.     if(m_prg_callback) {
  175.         m_prg_info.m_iter_total = N1*N2;
  176.         m_prg_info.m_iter_done = 0;
  177.         if(m_terminate = m_prg_callback(&m_prg_info)) {
  178.   return 0;
  179. }
  180.     }
  182.     bool bFreeGapLeft1  = m_esf_L1 && seg1 == m_Seq1;
  183.     bool bFreeGapRight1 = m_esf_R1 && m_Seq1 + m_SeqLen1 - len1 == seg1;
  184.     bool bFreeGapLeft2  = m_esf_L2 && seg2 == m_Seq2;
  185.     bool bFreeGapRight2 = m_esf_R2 && m_Seq2 + m_SeqLen2 - len2 == seg2;
  187.     TScore wgleft1   = bFreeGapLeft1? 0: m_Wg;
  188.     TScore wsleft1   = bFreeGapLeft1? 0: m_Ws;
  189.     TScore wg1 = m_Wg, ws1 = m_Ws;
  191.     // index calculation: [i,j] = i*n2 + j
  192.     vector<unsigned char> stl_bm (N1*N2);
  193.     unsigned char* backtrace_matrix = &stl_bm[0];
  195.     // first row
  196.     size_t k;
  197.     rowV[0] = wgleft1;
  198.     for (k = 1; k < N2; k++) {
  199.         rowV[k] = pV[k] + wsleft1;
  200.         rowF[k] = kInfMinus;
  201.         backtrace_matrix[k] = kMaskE | kMaskEc;
  202.     }
  203.     rowV[0] = 0;
  205.     if(m_prg_callback) {
  206.         m_prg_info.m_iter_done = k;
  207.         m_terminate = m_prg_callback(&m_prg_info);
  208.     }
  210.     // recurrences
  211.     TScore wgleft2   = bFreeGapLeft2? 0: m_Wg;
  212.     TScore wsleft2   = bFreeGapLeft2? 0: m_Ws;
  213.     TScore V  = rowV[N2 - 1];
  214.     TScore V0 = wgleft2;
  215.     TScore E, G, n0;
  216.     unsigned char tracer;
  218.     size_t i, j;
  219.     for(i = 1;  i < N1 && !m_terminate;  ++i) {
  220.         
  221.         V = V0 += wsleft2;
  222.         E = kInfMinus;
  223.         backtrace_matrix[k++] = kMaskFc;
  224.         unsigned char ci = seq1[i];
  226.         if(i == N1 - 1 && bFreeGapRight1) {
  227.                 wg1 = ws1 = 0;
  228.         }
  230.         TScore wg2 = m_Wg, ws2 = m_Ws;
  232.         for (j = 1; j < N2; ++j, ++k) {
  234.             G = pV[j] + sm[ci][(unsigned char)seq2[j]];
  235.             pV[j] = V;
  237.             n0 = V + wg1;
  238.             if(E >= n0) {
  239.                 E += ws1;      // continue the gap
  240.                 tracer = kMaskEc;
  241.             }
  242.             else {
  243.                 E = n0 + ws1;  // open a new gap
  244.                 tracer = 0;
  245.             }
  247.             if(j == N2 - 1 && bFreeGapRight2) {
  248.                 wg2 = ws2 = 0;
  249.             }
  250.             n0 = rowV[j] + wg2;
  251.             if(rowF[j] >= n0) {
  252.                 rowF[j] += ws2;
  253.                 tracer |= kMaskFc;
  254.             }
  255.             else {
  256.                 rowF[j] = n0 + ws2;
  257.             }
  259.             if (E >= rowF[j]) {
  260.                 if(E >= G) {
  261.                     V = E;
  262.                     tracer |= kMaskE;
  263.                 }
  264.                 else {
  265.                     V = G;
  266.                     tracer |= kMaskD;
  267.                 }
  268.             } else {
  269.                 if(rowF[j] >= G) {
  270.                     V = rowF[j];
  271.                 }
  272.                 else {
  273.                     V = G;
  274.                     tracer |= kMaskD;
  275.                 }
  276.             }
  277.             backtrace_matrix[k] = tracer;
  278.         }
  280.         pV[j] = V;
  282.         if(m_prg_callback) {
  283.             m_prg_info.m_iter_done = k;
  284.             if(m_terminate = m_prg_callback(&m_prg_info)) {
  285.                 break;
  286.             }
  287.         }
  288.        
  289.     }
  291.     if(!m_terminate) {
  292.         x_DoBackTrace(backtrace_matrix, N1, N2, transcript);
  293.     }
  295.     return V;
  296. }
  299. CNWAligner::TScore CNWAligner::Run()
  300. {
  301.     if(!m_Seq1 || !m_Seq2) {
  302.         NCBI_THROW(CAlgoAlignException, eNoData,
  303.                    "Sequence data not available for one or both sequences");
  304.     }
  306.     if(!x_CheckMemoryLimit()) {
  307.         NCBI_THROW(CAlgoAlignException, eMemoryLimit, "Out of space");
  308.     }
  310.     m_score = x_Run();
  312.     return m_score;
  313. }
  316. CNWAligner::TScore CNWAligner::x_Run()
  317. {
  318.     try {
  319.     
  320.         if(m_guides.size() == 0) {
  322.             m_score = x_Align(m_Seq1,m_SeqLen1, m_Seq2,m_SeqLen2,
  323.                               &m_Transcript);
  324.         }
  325.         else {
  326.             
  327.             // run the algorithm for every segment between hits
  328.             m_Transcript.clear();
  329.             size_t guides_dim = m_guides.size() / 4;
  330.             size_t q1 = m_SeqLen1, q0, s1 = m_SeqLen2, s0;
  331.             vector<ETranscriptSymbol> trans;
  332.             
  333.             // save original esf settings
  334.             bool esf_L1 = m_esf_L1, esf_R1 = m_esf_R1,
  335.                  esf_L2 = m_esf_L2, esf_R2 = m_esf_R2;
  336.             SetEndSpaceFree(false, esf_R1, false, esf_R2); // last region
  337.             for(size_t istart = 4*guides_dim, i = istart; i != 0; i -= 4) {
  339.                 q0 = m_guides[i - 3] + 1;
  340.                 s0 = m_guides[i - 1] + 1;
  341.                 size_t dim_query = q1 - q0, dim_subj = s1 - s0;
  342.                 x_Align(m_Seq1 + q0, dim_query, m_Seq2 + s0, dim_subj, &trans);
  343.                 copy(trans.begin(), trans.end(), back_inserter(m_Transcript));
  344.                 size_t dim_hit = m_guides[i - 3] - m_guides[i - 4] + 1;
  345.                 for(size_t k = 0; k < dim_hit; ++k) {
  346.                     m_Transcript.push_back(eTS_Match);
  347.                 }
  348.                 q1 = m_guides[i - 4];
  349.                 s1 = m_guides[i - 2];
  350.                 trans.clear();
  351.                 if(i == istart) {  // middle regions
  352.                     SetEndSpaceFree(false, false, false, false);
  353.                 }
  354.             }
  355.             SetEndSpaceFree(esf_L1, false, esf_L2, false); // first region
  356.             x_Align(m_Seq1, q1, m_Seq2, s1, &trans);
  357.             SetEndSpaceFree(esf_L1, esf_R1, esf_L2, esf_R2); // restore
  358.             copy(trans.begin(), trans.end(), back_inserter(m_Transcript));
  359.             m_score = x_ScoreByTranscript();
  360.         }
  361.     }
  362.     
  363.     catch( bad_alloc& ) {
  364.         
  365.         NCBI_THROW(CAlgoAlignException, eMemoryLimit, "Memory limit exceeded");
  366.     }
  367.     
  368.     return m_score;
  369. }
  372. // perform backtrace step;
  373. void CNWAligner::x_DoBackTrace(const unsigned char* backtrace,
  374.                                size_t N1, size_t N2,
  375.                                vector<ETranscriptSymbol>* transcript)
  376. {
  377.     transcript->clear();
  378.     transcript->reserve(N1 + N2);
  380.     size_t k = N1*N2 - 1;
  381.     while (k != 0) {
  382.         unsigned char Key = backtrace[k];
  383.         if (Key & kMaskD) {
  384.             transcript->push_back(eTS_Match);
  385.             k -= N2 + 1;
  386.         }
  387.         else if (Key & kMaskE) {
  388.             transcript->push_back(eTS_Insert); --k;
  389.             while(k > 0 && (Key & kMaskEc)) {
  390.                 transcript->push_back(eTS_Insert);
  391.                 Key = backtrace[k--];
  392.             }
  393.         }
  394.         else {
  395.             transcript->push_back(eTS_Delete);
  396.             k -= N2;
  397.             while(k > 0 && (Key & kMaskFc)) {
  398.                 transcript->push_back(eTS_Delete);
  399.                 Key = backtrace[k];
  400.                 k -= N2;
  401.             }
  402.         }
  403.     }
  404. }
  407. void  CNWAligner::SetPattern(const vector<size_t>& guides)
  408. {
  409.     size_t dim = guides.size();
  410.     const char* err = 0;
  411.     if(dim % 4 == 0) {
  412.         for(size_t i = 0; i < dim; i += 4) {
  414.             if( guides[i] > guides[i+1] || guides[i+2] > guides[i+3] ) {
  415.                 err = "Pattern hits must be specified in plus strand";
  416. break;
  417.             }
  419.             if(i > 4) {
  420.                 if(guides[i] <= guides[i-3] || guides[i+2] <= guides[i-2]){
  421.                     err = "Pattern hits coordinates must be sorted";
  422.                     break;
  423.                 }
  424.             }
  426.             size_t dim1 = guides[i + 1] - guides[i];
  427.             size_t dim2 = guides[i + 3] - guides[i + 2];
  428.             if( dim1 != dim2) {
  429.                 err = "Pattern hits must have equal length on both sequences";
  430.                 break;
  431.             }
  433.             if(guides[i+1] >= m_SeqLen1 || guides[i+3] >= m_SeqLen2) {
  434.                 err = "One or several pattern hits are out of range";
  435.                 break;
  436.             }
  437.         }
  438.     }
  439.     else {
  440.         err = "Pattern must have a dimension multiple of four";
  441.     }
  443.     if(err) {
  444.         NCBI_THROW(CAlgoAlignException, eBadParameter, err);
  445.     }
  446.     else {
  447.         m_guides = guides;
  448.     }
  449. }
  452. void CNWAligner::GetEndSpaceFree(bool* L1, bool* R1, bool* L2, bool* R2) const
  453. {
  454.     if(L1) *L1 = m_esf_L1;
  455.     if(R1) *R1 = m_esf_R1;
  456.     if(L2) *L2 = m_esf_L2;
  457.     if(R2) *R2 = m_esf_R2;
  458. }
  461. // Return transcript as a readable string
  462. string CNWAligner::GetTranscriptString(void) const
  463. {
  464.     const size_t dim = m_Transcript.size();   
  465.     string s;
  466.     s.resize(dim);
  467.     size_t i1 = 0, i2 = 0, i = 0;
  469.     for (int k = dim - 1; k >= 0;  --k) {
  471.         ETranscriptSymbol c0 = m_Transcript[k];
  472.         char c = 0;
  473.         switch ( c0 ) {
  475.             case eTS_Match: 
  476.             case eTS_Replace: {
  477.                 c = (m_Seq1[i1++] == m_Seq2[i2++])? 'M': 'R';
  478.             }
  479.             break;
  481.             case eTS_Insert: {
  482.                 c = 'I';
  483.                 ++i2;
  484.             }
  485.             break;
  487.             case eTS_Delete: {
  488.                 c = 'D';
  489.                 ++i1;
  490.             }
  491.             break;
  493.             case eTS_Intron: {
  494.                 c = '+';
  495.                 ++i2;
  496.             }
  497.             break;
  499.     default: {
  500.       NCBI_THROW(CAlgoAlignException, eInternal,
  501.                          "Unexpected transcript symbol");
  502.     }
  503.   
  504.         }
  506.         s[i++] = c;
  507.     }
  508.     if(i < s.size()) {
  509.         s.resize(i + 1);
  510.     }
  511.     return s;
  512. }
  515. void CNWAligner::SetProgressCallback ( FProgressCallback prg_callback,
  516.        void* data )
  517. {
  518.     m_prg_callback = prg_callback;
  519.     m_prg_info.m_data = data;
  520. }
  522.  
  523. void CNWAligner::SetScoreMatrix(const SNCBIPackedScoreMatrix* psm)
  524. {
  525.     if(psm) {
  526.         m_abc = psm->symbols;
  527. NCBISM_Unpack(psm, &m_ScoreMatrix);
  528.     }
  529.     else { // assume IUPACna
  530.         m_abc = g_nwaligner_nucleotides;
  531. memset(&m_ScoreMatrix.s, m_Wms, sizeof m_ScoreMatrix.s);
  532. // set diag to Wm except ambiguity symbols
  533.         unsigned char c1, c2;
  534. const size_t kNuclSize = strlen(g_nwaligner_nucleotides);
  535. for(size_t i = 0; i < kNuclSize; ++i) {
  536.     c1 = g_nwaligner_nucleotides[i];
  537.     for(size_t j = 0; j < i; ++j) {
  538.         c2 = g_nwaligner_nucleotides[j];
  539. m_ScoreMatrix.s[c1][c2] = m_ScoreMatrix.s[c2][c1] = m_Wms;
  540.     }
  541.     m_ScoreMatrix.s[c1][c1] = i < 4? m_Wm: m_Wms;
  542. }
  543.     }
  544. }
  547. // Check that all characters in sequence are valid for 
  548. // the current sequence type.
  549. // Return an index to the first invalid character
  550. // or len if all characters are valid.
  551. size_t CNWAligner::x_CheckSequence(const char* seq, size_t len) const
  552. {
  553.     char Flags [256];
  554.     memset(Flags, 0, sizeof Flags);
  555.     const size_t abc_size = strlen(m_abc);
  556.     
  557.     size_t k;
  558.     for(k = 0; k < abc_size; ++k) {
  559.         Flags[unsigned(m_abc[k])] = 1;
  560.     }
  562.     for(k = 0; k < len; ++k) {
  563.         if(Flags[unsigned(seq[k])] == 0)
  564.             break;
  565.     }
  567.     return k;
  568. }
  571. CNWAligner::TScore CNWAligner::GetScore() const 
  572. {
  573.   if(m_Transcript.size()) {
  574.       return m_score;
  575.   }
  576.   else {
  577.     NCBI_THROW(CAlgoAlignException, eNoData,
  578.                "Sequences not aligned yet");
  579.   }
  580. }
  583. bool CNWAligner::x_CheckMemoryLimit()
  584. {
  585.     const size_t elem_size = GetElemSize();
  586.     const size_t gdim = m_guides.size();
  587.     const size_t max_mem = kMax_UInt / 2;
  588.     if(gdim) {
  589.         size_t dim1 = m_guides[0], dim2 = m_guides[2];
  590.         double mem = double(dim1)*dim2*elem_size;
  591.         if(mem >= max_mem) {
  592.             return false;
  593.         }
  594.         for(size_t i = 4; i < gdim; i += 4) {
  595.             dim1 = m_guides[i] - m_guides[i-3] + 1;
  596.             dim2 = m_guides[i + 2] - m_guides[i-1] + 1;
  597.             mem = double(dim1)*dim2*elem_size;
  598.             if(mem >= max_mem) {
  599.                 return false;
  600.             }
  601.         }
  602.         dim1 = m_SeqLen1 - m_guides[gdim-3];
  603.         dim2 = m_SeqLen2 - m_guides[gdim-1];
  604.         mem = double(dim1)*dim2*elem_size;
  605.         if(mem >= max_mem) {
  606.             return false;
  607.         }
  609.         return true;
  610.     }
  611.     else {
  612.         double mem = double(m_SeqLen1 + 1)*(m_SeqLen2 + 1)*elem_size;
  613.         return mem < max_mem;
  614.     }
  615. }
  618. CNWAligner::TScore CNWAligner::x_ScoreByTranscript() const
  619. {
  620.     const size_t dim = m_Transcript.size();
  621.     if(dim == 0) return 0;
  623.     vector<ETranscriptSymbol> transcript (dim);
  624.     for(size_t i = 0; i < dim; ++i) {
  625.         transcript[i] = m_Transcript[dim - i - 1];
  626.     }
  628.     TScore score = 0;
  630.     const char* p1 = m_Seq1;
  631.     const char* p2 = m_Seq2;
  633.     int state1;   // 0 = normal, 1 = gap, 2 = intron
  634.     int state2;   // 0 = normal, 1 = gap
  636.     switch(transcript[0]) {
  637.     case eTS_Match:    state1 = state2 = 0; break;
  638.     case eTS_Insert:   state1 = 1; state2 = 0; score += m_Wg; break;
  639.     case eTS_Delete:   state1 = 0; state2 = 1; score += m_Wg; break;
  640.     default: {
  641.         NCBI_THROW(CAlgoAlignException, eInternal,
  642.                    "Invalid transcript symbol");
  643.         }
  644.     }
  646.     const TNCBIScore (*sm) [NCBI_FSM_DIM] = m_ScoreMatrix.s;
  648.     for(size_t i = 0; i < dim; ++i) {
  650.         unsigned char c1 = m_Seq1? *p1: 'N';
  651.         unsigned char c2 = m_Seq2? *p2: 'N';
  652.         switch(transcript[i]) {
  654.         case eTS_Match: {
  655.             state1 = state2 = 0;
  656.             score += sm[c1][c2];
  657.             ++p1; ++p2;
  658.         }
  659.         break;
  661.         case eTS_Insert: {
  662.             if(state1 != 1) score += m_Wg;
  663.             state1 = 1; state2 = 0;
  664.             score += m_Ws;
  665.             ++p2;
  666.         }
  667.         break;
  669.         case eTS_Delete: {
  670.             if(state2 != 1) score += m_Wg;
  671.             state1 = 0; state2 = 1;
  672.             score += m_Ws;
  673.             ++p1;
  674.         }
  675.         break;
  676.         
  677.         default: {
  678.         NCBI_THROW(CAlgoAlignException, eInternal,
  679.                    "Invalid transcript symbol");
  680.         }
  681.         }
  682.     }
  684.     if(m_esf_L1) {
  685.         size_t g = 0;
  686.         for(size_t i = 0; i < dim; ++i) {
  687.             if(transcript[i] == eTS_Insert) ++g; else break;
  688.         }
  689.         if(g > 0) {
  690.             score -= (m_Wg + g*m_Ws);
  691.         }
  692.     }
  694.     if(m_esf_L2) {
  695.         size_t g = 0;
  696.         for(size_t i = 0; i < dim; ++i) {
  697.             if(transcript[i] == eTS_Delete) ++g; else break;
  698.         }
  699.         if(g > 0) {
  700.             score -= (m_Wg + g*m_Ws);
  701.         }
  702.     }
  704.     if(m_esf_R1) {
  705.         size_t g = 0;
  706.         for(int i = dim - 1; i >= 0; --i) {
  707.             if(transcript[i] == eTS_Insert) ++g; else break;
  708.         }
  709.         if(g > 0) {
  710.             score -= (m_Wg + g*m_Ws);
  711.         }
  712.     }
  714.     if(m_esf_R2) {
  715.         size_t g = 0;
  716.         for(int i = dim - 1; i >= 0; --i) {
  717.             if(transcript[i] == eTS_Delete) ++g; else break;
  718.         }
  719.         if(g > 0) {
  720.             score -= (m_Wg + g*m_Ws);
  721.         }
  722.     }
  724.     return score;
  725. }
  728. size_t CNWAligner::GetLeftSeg(size_t* q0, size_t* q1,
  729.                               size_t* s0, size_t* s1,
  730.                               size_t min_size) const
  731. {
  732.     size_t trdim = m_Transcript.size();
  733.     size_t cur = 0, maxseg = 0;
  734.     const char* p1 = m_Seq1;
  735.     const char* p2 = m_Seq2;
  736.     size_t i0 = 0, j0 = 0, imax = i0, jmax = j0;
  738.     for(int k = trdim - 1; k >= 0; --k) {
  740.         switch(m_Transcript[k]) {
  742.             case eTS_Insert: {
  743.                 ++p2;
  744.                 if(cur > maxseg) {
  745.                     maxseg = cur;
  746.                     imax = i0;
  747.                     jmax = j0;
  748.                     if(maxseg >= min_size) goto ret_point;
  749.                 }
  750.                 cur = 0;
  751.             }
  752.             break;
  754.             case eTS_Delete: {
  755.             ++p1;
  756.             if(cur > maxseg) {
  757.                 maxseg = cur;
  758.                 imax = i0;
  759.                 jmax = j0;
  760.                 if(maxseg >= min_size) goto ret_point;
  761.             }
  762.             cur = 0;
  763.             }
  764.             break;
  766.             case eTS_Match:
  767.             case eTS_Replace: {
  768.                 if(*p1 == *p2) {
  769.                     if(cur == 0) {
  770.                         i0 = p1 - m_Seq1;
  771.                         j0 = p2 - m_Seq2;
  772.                     }
  773.                     ++cur;
  774.                 }
  775.                 else {
  776.                     if(cur > maxseg) {
  777.                         maxseg = cur;
  778.                         imax = i0;
  779.                         jmax = j0;
  780.                         if(maxseg >= min_size) goto ret_point;
  781.                     }
  782.                     cur = 0;
  783.                 }
  784.                 ++p1;
  785.                 ++p2;
  786.             }
  787.             break;
  789.             default: {
  790.                 NCBI_THROW(CAlgoAlignException, eInternal,
  791.                            "Invalid transcript symbol");
  792.             }
  793.         }
  794.     }
  796.     if(cur > maxseg) {
  797.       maxseg = cur;
  798.       imax = i0;
  799.       jmax = j0;
  800.     }
  802.  ret_point:
  804.     *q0 = imax; *s0 = jmax;
  805.     *q1 = *q0 + maxseg - 1;
  806.     *s1 = *s0 + maxseg - 1;
  808.     return maxseg;
  809. }
  812. size_t CNWAligner::GetRightSeg(size_t* q0, size_t* q1,
  813.                                size_t* s0, size_t* s1,
  814.                                size_t min_size) const
  815. {
  816.     size_t trdim = m_Transcript.size();
  817.     size_t cur = 0, maxseg = 0;
  818.     const char* seq1_end = m_Seq1 + m_SeqLen1;
  819.     const char* seq2_end = m_Seq2 + m_SeqLen2;
  820.     const char* p1 = seq1_end - 1;
  821.     const char* p2 = seq2_end - 1;
  822.     size_t i0 = m_SeqLen1 - 1, j0 = m_SeqLen2 - 1,
  823.            imax = i0, jmax = j0;
  825.     for(size_t k = 0; k < trdim; ++k) {
  827.         switch(m_Transcript[k]) {
  829.             case eTS_Insert: {
  830.                 --p2;
  831.                 if(cur > maxseg) {
  832.                     maxseg = cur;
  833.                     imax = i0;
  834.                     jmax = j0;
  835.                     if(maxseg >= min_size) goto ret_point;
  836.                 }
  837.                 cur = 0;
  838.             }
  839.             break;
  841.             case eTS_Delete: {
  842.                 --p1;
  843.                 if(cur > maxseg) {
  844.                     maxseg = cur;
  845.                     imax = i0;
  846.                     jmax = j0;
  847.                     if(maxseg >= min_size) goto ret_point;
  848. }
  849.                 cur = 0;
  850.             }
  851.             break;
  853.             case eTS_Match:
  854.             case eTS_Replace: {
  855.                 if(*p1 == *p2) {
  856.                     if(cur == 0) {
  857.                         i0 = p1 - m_Seq1;
  858.                         j0 = p2 - m_Seq2;
  859.                     }
  860.                     ++cur;
  861.                 }
  862.                 else {
  863.                     if(cur > maxseg) {
  864.                         maxseg = cur;
  865.                         imax = i0;
  866.                         jmax = j0;
  867.                         if(maxseg >= min_size) goto ret_point;
  868.                     }
  869.                     cur = 0;
  870.                 }
  871.                 --p1;
  872.                 --p2;
  873.             }
  874.             break;
  876.             default: {
  877.                 NCBI_THROW(CAlgoAlignException, eInternal,
  878.                            "Invalid transcript symbol");
  879.             }
  880.         }
  881.     }
  883.     if(cur > maxseg) {
  884.       maxseg = cur;
  885.       imax = i0;
  886.       jmax = j0;
  887.     }
  889.  ret_point:
  891.     *q1 = imax; *s1 = jmax;
  892.     *q0 = imax - maxseg + 1;
  893.     *s0 = jmax - maxseg + 1;
  895.     return maxseg;
  896. }
  899. size_t CNWAligner::GetLongestSeg(size_t* q0, size_t* q1,
  900.                                  size_t* s0, size_t* s1) const
  901. {
  902.     size_t trdim = m_Transcript.size();
  903.     size_t cur = 0, maxseg = 0;
  904.     const char* p1 = m_Seq1;
  905.     const char* p2 = m_Seq2;
  906.     size_t i0 = 0, j0 = 0, imax = i0, jmax = j0;
  908.     for(int k = trdim-1; k >= 0; --k) {
  910.         switch(m_Transcript[k]) {
  912.             case eTS_Insert: {
  913.                 ++p2;
  914.                 if(cur > maxseg) {
  915.                     maxseg = cur;
  916.                     imax = i0;
  917.                     jmax = j0;
  918.                 }
  919.                 cur = 0;
  920.             }
  921.             break;
  923.             case eTS_Delete: {
  924.             ++p1;
  925.             if(cur > maxseg) {
  926.                 maxseg = cur;
  927.                 imax = i0;
  928.                 jmax = j0;
  929.             }
  930.             cur = 0;
  931.             }
  932.             break;
  934.             case eTS_Match:
  935.             case eTS_Replace: {
  936.                 if(*p1 == *p2) {
  937.                     if(cur == 0) {
  938.                         i0 = p1 - m_Seq1;
  939.                         j0 = p2 - m_Seq2;
  940.                     }
  941.                     ++cur;
  942.                 }
  943.                 else {
  944.     if(cur > maxseg) {
  945.                         maxseg = cur;
  946.                         imax = i0;
  947.                         jmax = j0;
  948.                     }
  949.                     cur = 0;
  950.                 }
  951.                 ++p1;
  952.                 ++p2;
  953.             }
  954.             break;
  956.             default: {
  957.                 NCBI_THROW(CAlgoAlignException, eInternal,
  958.                            "Invalid transcript symbol");
  959.             }
  960.         }
  961.     }
  962.     
  963.     if(cur > maxseg) {
  964.         maxseg = cur;
  965.         imax = i0;
  966.         jmax = j0;
  967.     }
  968.     
  969.     *q0 = imax; *s0 = jmax;
  970.     *q1 = *q0 + maxseg - 1;
  971.     *s1 = *s0 + maxseg - 1;
  973.     return maxseg;
  974. }
  977. END_NCBI_SCOPE
  980. /*
  981.  * ===========================================================================
  982.  * $Log: nw_aligner.cpp,v $
  983.  * Revision 1000.3  2004/06/01 18:04:52  gouriano
  984.  * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.48
  985.  *
  986.  * Revision 1.48  2004/05/21 21:41:02  gorelenk
  987.  * Added PCH ncbi_pch.hpp
  988.  *
  989.  * Revision 1.47  2004/05/19 13:37:47  kapustin
  990.  * Remove test dumping code
  991.  *
  992.  * Revision 1.46  2004/05/18 21:43:40  kapustin
  993.  * Code cleanup
  994.  *
  995.  * Revision 1.45  2004/05/17 14:50:56  kapustin
  996.  * Add/remove/rearrange some includes and object declarations
  997.  *
  998.  * Revision 1.44  2004/04/23 14:39:47  kapustin
  999.  * Add Splign library and other changes
  1000.  *
  1001.  * Revision 1.42  2003/12/29 13:03:48  kapustin
  1002.  * Return string from GetTranscriptString().
  1003.  *
  1004.  * Revision 1.41  2003/10/31 19:40:13  kapustin
  1005.  * Get rid of some WS and GCC complains
  1006.  *
  1007.  * Revision 1.40  2003/10/27 20:47:35  kapustin
  1008.  * Minor code cleanup
  1009.  *
  1010.  * Revision 1.39  2003/10/14 18:44:46  kapustin
  1011.  * Adjust esf flags during pattern-guided alignment
  1012.  *
  1013.  * Revision 1.38  2003/09/30 19:50:04  kapustin
  1014.  * Make use of standard score matrix interface
  1015.  *
  1016.  * Revision 1.37  2003/09/26 14:43:18  kapustin
  1017.  * Remove exception specifications
  1018.  *
  1019.  * Revision 1.36  2003/09/15 20:49:11  kapustin
  1020.  * Clear the transcript when setting new sequences
  1021.  *
  1022.  * Revision 1.35  2003/09/04 16:07:38  kapustin
  1023.  * Use STL vectors for exception-safe dynamic arrays and matrices
  1024.  *
  1025.  * Revision 1.34  2003/09/03 17:28:40  kapustin
  1026.  * Clean the list of includes
  1027.  *
  1028.  * Revision 1.33  2003/09/02 22:38:18  kapustin
  1029.  * Move format functionality out of the class
  1030.  *
  1031.  * Revision 1.32  2003/07/09 15:11:22  kapustin
  1032.  * Update plain text output with seq ids and use inequalities instead of
  1033.  * binary operation to verify intron boundaries
  1034.  *
  1035.  * Revision 1.31  2003/06/26 20:39:39  kapustin
  1036.  * Rename formal parameters in SetEndSpaceFree() to avoid conflict with macro
  1037.  * under some configurations
  1038.  *
  1039.  * Revision 1.30  2003/06/17 17:20:44  kapustin
  1040.  * CNWAlignerException -> CAlgoAlignException
  1041.  *
  1042.  * Revision 1.29  2003/06/17 16:06:13  kapustin
  1043.  * Detect all variety of splices in Seq-Align formatter (restored from 1.27)
  1044.  *
  1045.  * Revision 1.28  2003/06/17 14:51:04  dicuccio
  1046.  * Fixed after algo/ rearragnement
  1047.  *
  1048.  * Revision 1.26  2003/06/02 14:04:49  kapustin
  1049.  * Progress indication-related updates
  1050.  *
  1051.  * Revision 1.25  2003/05/23 18:27:02  kapustin
  1052.  * Use weak comparisons in core recurrences. Adjust for new
  1053.  * transcript identifiers.
  1054.  *
  1055.  * Revision 1.24  2003/04/17 14:44:40  kapustin
  1056.  * A few changes to eliminate gcc warnings
  1057.  *
  1058.  * Revision 1.23  2003/04/14 19:00:55  kapustin
  1059.  * Add guide creation facility.  x_Run() -> x_Align()
  1060.  *
  1061.  * Revision 1.22  2003/04/10 19:15:16  kapustin
  1062.  * Introduce guiding hits approach
  1063.  *
  1064.  * Revision 1.21  2003/04/02 20:52:55  kapustin
  1065.  * Make FormatAsText virtual. Pass output string as a parameter.
  1066.  *
  1067.  * Revision 1.20  2003/03/31 15:32:05  kapustin
  1068.  * Calculate score independently from transcript
  1069.  *
  1070.  * Revision 1.19  2003/03/19 16:36:39  kapustin
  1071.  * Reverse memory limit condition
  1072.  *
  1073.  * Revision 1.18  2003/03/18 15:15:51  kapustin
  1074.  * Implement virtual memory checking function. Allow separate free end
  1075.  * gap specification
  1076.  *
  1077.  * Revision 1.17  2003/03/17 15:31:46  kapustin
  1078.  * Forced conversion to double in memory limit checking to avoid integer
  1079.  * overflow
  1080.  *
  1081.  * Revision 1.16  2003/03/14 19:18:50  kapustin
  1082.  * Add memory limit checking. Fix incorrect seq index references when
  1083.  * reporting about incorrect input seq character. Support all characters
  1084.  * within IUPACna coding
  1085.  *
  1086.  * Revision 1.15  2003/03/07 13:51:11  kapustin
  1087.  * Use zero-based indices to specify seq coordinates in ASN
  1088.  *
  1089.  * Revision 1.14  2003/03/05 20:13:48  kapustin
  1090.  * Simplify FormatAsText(). Fix FormatAsSeqAlign(). Convert sequence
  1091.  * alphabets to capitals
  1092.  *
  1093.  * Revision 1.13  2003/02/11 16:06:54  kapustin
  1094.  * Add end-space free alignment support
  1095.  *
  1096.  * Revision 1.12  2003/02/04 23:04:50  kapustin
  1097.  * Add progress callback support
  1098.  *
  1099.  * Revision 1.11  2003/01/31 02:55:48  lavr
  1100.  * User proper argument for ::time()
  1101.  *
  1102.  * Revision 1.10  2003/01/30 20:32:36  kapustin
  1103.  * Add EstiamteRunningTime()
  1104.  *
  1105.  * Revision 1.9  2003/01/28 12:39:03  kapustin
  1106.  * Implement ASN.1 and SeqAlign output
  1107.  *
  1108.  * Revision 1.8  2003/01/24 16:48:50  kapustin
  1109.  * Support more output formats - type 2 and gapped FastA
  1110.  *
  1111.  * Revision 1.7  2003/01/21 16:34:21  kapustin
  1112.  * Get rid of reverse() and reverse_copy() not supported by MSVC and/or ICC
  1113.  *
  1114.  * Revision 1.6  2003/01/21 12:41:37  kapustin
  1115.  * Use class neg infinity constant to specify least possible score
  1116.  *
  1117.  * Revision 1.5  2003/01/08 15:42:59  kapustin
  1118.  * Fix initialization for the first column of the backtrace matrix
  1119.  *
  1120.  * Revision 1.4  2003/01/06 18:04:10  kapustin
  1121.  * CNWAligner: add sequence size verification
  1122.  *
  1123.  * Revision 1.3  2002/12/31 13:53:13  kapustin
  1124.  * CNWAligner::Format() -- use CNcbiOstrstreamToString
  1125.  *
  1126.  * Revision 1.2  2002/12/12 17:59:28  kapustin
  1127.  * Enable spliced alignments
  1128.  *
  1129.  * Revision 1.1  2002/12/06 17:41:21  ivanov
  1130.  * Initial revision
  1131.  *
  1132.  * ===========================================================================
  1133.  */