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

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

生物技术

开发平台：

C/C++

nw_formatter.cpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: nw_formatter.cpp,v $
  4.  * PRODUCTION Revision 1000.2  2004/06/01 18:04:54  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.10
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /* $Id: nw_formatter.cpp,v 1000.2 2004/06/01 18:04:54 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.  * Author:  Yuri Kapustin
  36.  *
  37.  * ===========================================================================
  38.  *
  39.  */
  41. #include <ncbi_pch.hpp>
  42. #include <algo/align/nw_formatter.hpp>
  43. #include <algo/align/nw_aligner.hpp>
  44. #include <algo/align/align_exception.hpp>
  46. #include <objects/seqalign/Score.hpp>
  47. #include <objects/general/Object_id.hpp>
  48. #include <objects/seqalign/Dense_seg.hpp>
  49. #include <objects/seqloc/Seq_id.hpp>
  50. #include <objects/seqalign/Seq_align.hpp>
  51. #include <serial/objostrasn.hpp>
  52. #include <serial/serial.hpp>
  54. #include <iterator>
  56. BEGIN_NCBI_SCOPE
  57. USING_SCOPE(objects);
  60. CNWFormatter::CNWFormatter (const CNWAligner& aligner):
  61.     m_aligner(&aligner)
  62. {
  63. }
  66. void CNWFormatter::AsSeqAlign(CSeq_align* seqalign) const
  67. {
  68.     if(seqalign == 0) {
  69.         NCBI_THROW(CAlgoAlignException,
  70.                    eBadParameter,
  71.                    "Invalid address specified");
  72.     }
  73.     
  74.     const vector<CNWAligner::ETranscriptSymbol>& transcript = 
  75.         *(m_aligner->GetTranscript());
  77.     if(transcript.size() == 0) {
  78.         NCBI_THROW(CAlgoAlignException,
  79.                    eNoData,
  80.                    "Zero size transcript (forgot to run the aligner?)");
  81.     }
  83.     seqalign->Reset();
  85.     // the alignment is pairwise
  86.     seqalign->SetDim(2);
  88.     // this is a global alignment
  89.     seqalign->SetType(CSeq_align::eType_global);
  91.     // seq-ids
  92.     CRef< CSeq_id > id1 ( new CSeq_id );
  93.     CRef< CObject_id > local_oid1 (new CObject_id);
  94.     local_oid1->SetStr(m_Seq1Id);
  95.     id1->SetLocal(*local_oid1);
  97.     CRef< CSeq_id > id2 ( new CSeq_id );
  98.     CRef< CObject_id > local_oid2 (new CObject_id);
  99.     local_oid2->SetStr(m_Seq2Id);
  100.     id2->SetLocal(*local_oid2);
  101.     
  102.     // the score was calculated during the main process
  103.     CRef< CScore > score (new CScore);
  104.     CRef< CObject_id > id (new CObject_id);
  105.     id->SetStr("score");
  106.     score->SetId(*id);
  107.     CRef< CScore::C_Value > val (new CScore::C_Value);
  108.     val->SetInt(m_aligner->GetScore());
  109.     score->SetValue(*val);
  110.     CSeq_align::TScore& scorelist = seqalign->SetScore();
  111.     scorelist.push_back(score);
  113.     // create segments and add them to this seq-align
  114.     CRef< CSeq_align::C_Segs > segs (new CSeq_align::C_Segs);
  115.     CDense_seg& ds = segs->SetDenseg();
  116.     ds.SetDim(2);
  117.     CDense_seg::TIds& ids = ds.SetIds();
  118.     ids.push_back( id1 );
  119.     ids.push_back( id2 );
  120.     CDense_seg::TStarts&  starts  = ds.SetStarts();
  121.     CDense_seg::TLens&    lens    = ds.SetLens();
  122.     CDense_seg::TStrands& strands = ds.SetStrands();
  123.     
  124.     // iterate through transcript
  125.     size_t seg_count = 0;
  126.     {{ 
  127.         const char * const S1 = m_aligner->GetSeq1();
  128.         const char * const S2 = m_aligner->GetSeq2();
  129.         const char *seq1 = S1, *seq2 = S2;
  130.         const char *start1 = seq1, *start2 = seq2;
  132.         vector<CNWAligner::ETranscriptSymbol>::const_reverse_iterator
  133.             ib = transcript.rbegin(),
  134.             ie = transcript.rend(),
  135.             ii;
  136.         
  137.         CNWAligner::ETranscriptSymbol ts = *ib;
  138.         bool intron = (ts == CNWAligner::eTS_Intron);
  139.         char seg_type0 = ((ts == CNWAligner::eTS_Insert || intron )? 1:
  140.                           (ts == CNWAligner::eTS_Delete)? 2: 0);
  141.         size_t seg_len = 0;
  143.         for (ii = ib;  ii != ie; ++ii) {
  144.             ts = *ii;
  145.             intron = (ts == CNWAligner::eTS_Intron);
  146.             char seg_type = ((ts == CNWAligner::eTS_Insert || intron )? 1:
  147.                              (ts == CNWAligner::eTS_Delete)? 2: 0);
  148.             if(seg_type0 != seg_type) {
  149.                 starts.push_back( (seg_type0 == 1)? -1: start1 - S1 );
  150.                 starts.push_back( (seg_type0 == 2)? -1: start2 - S2 );
  151.                 lens.push_back(seg_len);
  152.                 strands.push_back(eNa_strand_plus);
  153.                 strands.push_back(eNa_strand_plus);
  154.                 ++seg_count;
  155.                 start1 = seq1;
  156.                 start2 = seq2;
  157.                 seg_type0 = seg_type;
  158.                 seg_len = 1;
  159.             }
  160.             else {
  161.                 ++seg_len;
  162.             }
  164.             if(seg_type != 1) ++seq1;
  165.             if(seg_type != 2) ++seq2;
  166.         }
  167.         // the last one
  168.         starts.push_back( (seg_type0 == 1)? -1: start1 - S1 );
  169.         starts.push_back( (seg_type0 == 2)? -1: start2 - S2 );
  170.         lens.push_back(seg_len);
  171.         strands.push_back(eNa_strand_plus);
  172.         strands.push_back(eNa_strand_plus);
  173.         ++seg_count;
  174.     }}
  176.     ds.SetNumseg(seg_count);
  177.     ds.SetIds();
  178.     seqalign->SetSegs(*segs);
  179. }
  183. void CNWFormatter::AsText(string* output, ETextFormatType type,
  184.                           size_t line_width) const
  185. {
  186.     CNcbiOstrstream ss;
  188.     const vector<CNWAligner::ETranscriptSymbol>& transcript =
  189.         *(m_aligner->GetTranscript());
  191.     if(transcript.size() == 0) {
  192.         NCBI_THROW(CAlgoAlignException,
  193.                    eNoData,
  194.                    "Zero size transcript (forgot to run the aligner?)");
  195.     }
  197.     switch (type) {
  199.     case eFormatType1: {
  200.         if(m_Seq1Id.size() && m_Seq2Id.size()) {
  201.             ss << '>' << m_Seq1Id << 't' << m_Seq2Id << endl;
  202.         }
  203.         vector<char> v1, v2;
  204.         size_t aln_size = x_ApplyTranscript(&v1, &v2);
  205.         unsigned i1 = 0, i2 = 0;
  206.         for (size_t i = 0;  i < aln_size; ) {
  207.             ss << i << 't' << i1 << ':' << i2 << endl;
  208.             int i0 = i;
  209.             for (size_t jPos = 0;  i < aln_size  &&  jPos < line_width;
  210.                  ++i, ++jPos) {
  211.                 char c = v1[i0 + jPos];
  212.                 ss << c;
  213.                 if(c != '-'  &&  c != '+')
  214.                     i1++;
  215.             }
  216.             ss << endl;
  217.             
  218.             string marker_line(line_width, ' ');
  219.             i = i0;
  220.             for (size_t jPos = 0;  i < aln_size  &&  jPos < line_width;
  221.                  ++i, ++jPos) {
  222.                 char c1 = v1[i0 + jPos];
  223.                 char c  = v2[i0 + jPos];
  224.                 ss << c;
  225.                 if(c != '-' && c != '+')
  226.                     i2++;
  227.                 if(c != c1  &&  c != '-'  &&  c1 != '-'  &&  c1 != '+')
  228.                     marker_line[jPos] = '^';
  229.             }
  230.             ss << endl << marker_line << endl;
  231.         }
  232.     }
  233.     break;
  235.     case eFormatType2: {
  236.         if(m_Seq1Id.size() && m_Seq2Id.size()) {
  237.             ss << '>' << m_Seq1Id << 't' << m_Seq2Id << endl;
  238.         }
  239.         vector<char> v1, v2;
  240.         size_t aln_size = x_ApplyTranscript(&v1, &v2);
  241.         unsigned i1 = 0, i2 = 0;
  242.         for (size_t i = 0;  i < aln_size; ) {
  243.             ss << i << 't' << i1 << ':' << i2 << endl;
  244.             int i0 = i;
  245.             for (size_t jPos = 0;  i < aln_size  &&  jPos < line_width;
  246.                  ++i, ++jPos) {
  247.                 char c = v1[i0 + jPos];
  248.                 ss << c;
  249.                 if(c != '-'  &&  c != '+')
  250.                     i1++;
  251.             }
  252.             ss << endl;
  253.             
  254.             string line2 (line_width, ' ');
  255.             string line3 (line_width, ' ');
  256.             i = i0;
  257.             for (size_t jPos = 0;  i < aln_size  &&  jPos < line_width;
  258.                  ++i, ++jPos) {
  259.                 char c1 = v1[i0 + jPos];
  260.                 char c2  = v2[i0 + jPos];
  261.                 if(c2 != '-' && c2 != '+')
  262.                     i2++;
  263.                 if(c2 == c1)
  264.                     line2[jPos] = '|';
  265.                 line3[jPos] = c2;
  266.             }
  267.             ss << line2 << endl << line3 << endl << endl;
  268.         }
  269.     }
  270.     break;
  272.     case eFormatAsn: {
  273.         CSeq_align seq_align;
  274.         AsSeqAlign(&seq_align);
  275.         CObjectOStreamAsn asn_stream (ss);
  276.         asn_stream << seq_align;
  277.         asn_stream << Separator;
  278.     }
  279.     break;
  281.     case eFormatFastA: {
  282.         vector<char> v1, v2;
  283.         size_t aln_size = x_ApplyTranscript(&v1, &v2);
  284.         
  285.         ss << '>' << m_Seq1Id << endl;
  286.         const vector<char>* pv = &v1;
  287.         for(size_t i = 0; i < aln_size; ++i) {
  288.             for(size_t j = 0; j < line_width && i < aln_size; ++j, ++i) {
  289.                 ss << (*pv)[i];
  290.             }
  291.             ss << endl;
  292.         }
  294.         ss << '>' << m_Seq2Id << endl;
  295.         pv = &v2;
  296.         for(size_t i = 0; i < aln_size; ++i) {
  297.             for(size_t j = 0; j < line_width && i < aln_size; ++j, ++i) {
  298.                 ss << (*pv)[i];
  299.             }
  300.             ss << endl;
  301.         }
  302.     }
  303.     break;
  304.     
  305.     case eFormatExonTable:
  306.     case eFormatExonTableEx: {
  307.         ss.precision(3);
  308.         bool esfL1, esfR1, esfL2, esfR2;
  309.         m_aligner->GetEndSpaceFree(&esfL1, &esfR1, &esfL2, &esfR2);
  310.         const size_t len2 = m_aligner->GetSeqLen2();
  311.         const char* start1 = m_aligner->GetSeq1();
  312.         const char* start2 = m_aligner->GetSeq2();
  313.         const char* p1 = start1;
  314.         const char* p2 = start2;
  315.         int tr_idx_hi0 = transcript.size() - 1, tr_idx_hi = tr_idx_hi0;
  316.         int tr_idx_lo0 = 0, tr_idx_lo = tr_idx_lo0;
  318.         if(esfL1 && transcript[tr_idx_hi0] == CNWAligner::eTS_Insert) {
  319.             while(esfL1 && transcript[tr_idx_hi] == CNWAligner::eTS_Insert) {
  320.                 --tr_idx_hi;
  321.                 ++p2;
  322.             }
  323.         }
  325.         if(esfL2 && transcript[tr_idx_hi0] == CNWAligner::eTS_Delete) {
  326.             while(esfL2 && transcript[tr_idx_hi] == CNWAligner::eTS_Delete) {
  327.                 --tr_idx_hi;
  328.                 ++p1;
  329.             }
  330.         }
  332.         if(esfR1 && transcript[tr_idx_lo0] == CNWAligner::eTS_Insert) {
  333.             while(esfR1 && transcript[tr_idx_lo] == CNWAligner::eTS_Insert) {
  334.                 ++tr_idx_lo;
  335.             }
  336.         }
  338.         if(esfR2 && transcript[tr_idx_lo0] == CNWAligner::eTS_Delete) {
  339.             while(esfR2 && transcript[tr_idx_lo] == CNWAligner::eTS_Delete) {
  340.                 ++tr_idx_lo;
  341.             }
  342.         }
  344.         bool type_ex = type == eFormatExonTableEx;
  345.         vector<char> trans_ex (tr_idx_hi - tr_idx_lo + 1);
  347.         for(int tr_idx = tr_idx_hi; tr_idx >= tr_idx_lo;) {
  348.             const char* p1_beg = p1;
  349.             const char* p2_beg = p2;
  350.             size_t matches = 0, exon_aln_size = 0;
  352.             vector<char>::iterator ii_ex = trans_ex.begin();
  353.             while(tr_idx >= tr_idx_lo && 
  354.                   transcript[tr_idx] < CNWAligner::eTS_Intron) {
  355.                 
  356.                 bool noins = transcript[tr_idx] != CNWAligner::eTS_Insert;
  357.                 bool nodel = transcript[tr_idx] != CNWAligner::eTS_Delete;
  358.                 if(noins && nodel) {
  359.                     if(*p1++ == *p2++) {
  360.                         ++matches;
  361.                         if(type_ex) *ii_ex++ = 'M';
  362.                     }
  363.                     else {
  364.                         if(type_ex) *ii_ex++ = 'R';
  365.                     }
  366.                 } else if( noins ) {
  367.                     ++p1;
  368.                     if(type_ex) *ii_ex++ = 'D';
  369.                 } else {
  370.                     ++p2;
  371.                     if(type_ex) *ii_ex++ = 'I';
  372.                 }
  373.                 --tr_idx;
  374.                 ++exon_aln_size;
  375.             }
  377.             if(exon_aln_size > 0) {
  379.                 if(m_Seq1Id.size() && m_Seq2Id.size()) {
  380.                     ss << m_Seq1Id << 't' << m_Seq2Id << 't';
  381.                 }
  382. else {
  383.     ss << "-t-t";
  384. }
  385.                 float identity = float(matches) / exon_aln_size;
  386.                 ss << identity << 't' << exon_aln_size << 't';
  387.                 size_t beg1  = p1_beg - start1, end1 = p1 - start1 - 1;
  388.                 size_t beg2  = p2_beg - start2, end2 = p2 - start2 - 1;
  389. if(beg1 <= end1) {
  390.     ss << beg1 << 't' << end1 << 't';
  391. }
  392. else {
  393.     ss << "-t-t";
  394. }
  395. ss << beg2 << 't' << end2 << 't';
  396.                 char c1 = (p2_beg >= start2 + 2)? *(p2_beg - 2): ' ';
  397.                 char c2 = (p2_beg >= start2 + 1)? *(p2_beg - 1): ' ';
  398.                 char c3 = (p2 < start2 + len2)? *(p2): ' ';
  399.                 char c4 = (p2 < start2 + len2 - 1)? *(p2+1): ' ';
  400.                 ss << c1 << c2 << "<exon>" << c3 << c4;
  402.                 if(type == eFormatExonTableEx) {
  403.                     ss << 't';
  404.                     copy(trans_ex.begin(), ii_ex, ostream_iterator<char>(ss));
  405.                 }
  406.                 ss << endl;
  407.             }
  409.             // find next exon
  410.             while(tr_idx >= tr_idx_lo &&
  411.                   (transcript[tr_idx] == CNWAligner::eTS_Intron)) {
  413.                 --tr_idx;
  414.                 ++p2;
  415.             }
  416.         }
  417.     }
  418.     break;
  420.     default: {
  421.         NCBI_THROW(CAlgoAlignException, eBadParameter,
  422.                    "Incorrect format specified");
  423.     }
  425.     }
  427.     *output = CNcbiOstrstreamToString(ss);
  428. }
  432. // Transform source sequences according to the transcript.
  433. // Write the results to v1 and v2 leaving source sequences intact.
  434. // Return alignment size.
  435. size_t CNWFormatter::x_ApplyTranscript(vector<char>* pv1, vector<char>* pv2)
  436.     const
  437. {
  438.     const vector<CNWAligner::ETranscriptSymbol>& transcript = 
  439.         *(m_aligner->GetTranscript());
  441.     vector<char>& v1 = *pv1;
  442.     vector<char>& v2 = *pv2;
  444.     vector<CNWAligner::ETranscriptSymbol>::const_reverse_iterator
  445.         ib = transcript.rbegin(),
  446.         ie = transcript.rend(),
  447.         ii;
  449.     const char* iv1 = m_aligner->GetSeq1();
  450.     const char* iv2 = m_aligner->GetSeq2();
  451.     v1.clear();
  452.     v2.clear();
  454.     for (ii = ib;  ii != ie;  ii++) {
  455.         CNWAligner::ETranscriptSymbol ts = *ii;
  456.         char c1, c2;
  457.         switch ( ts ) {
  458.         case CNWAligner::eTS_Insert:
  459.             c1 = '-';
  460.             c2 = *iv2++;
  461.             break;
  462.         case CNWAligner::eTS_Delete:
  463.             c2 = '-';
  464.             c1 = *iv1++;
  465.             break;
  466.         case CNWAligner::eTS_Match:
  467.         case CNWAligner::eTS_Replace:
  468.             c1 = *iv1++;
  469.             c2 = *iv2++;
  470.             break;
  471.         case CNWAligner::eTS_Intron:
  472.             c1 = '+';
  473.             c2 = *iv2++;
  474.             break;
  475.         default:
  476.             c1 = c2 = '?';
  477.             break;
  478.         }
  479.         v1.push_back(c1);
  480.         v2.push_back(c2);
  481.     }
  483.     return v1.size();
  484. }
  487. END_NCBI_SCOPE
  490. /*
  491.  * ===========================================================================
  492.  * $Log: nw_formatter.cpp,v $
  493.  * Revision 1000.2  2004/06/01 18:04:54  gouriano
  494.  * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.10
  495.  *
  496.  * Revision 1.10  2004/05/21 21:41:02  gorelenk
  497.  * Added PCH ncbi_pch.hpp
  498.  *
  499.  * Revision 1.9  2004/05/18 21:43:40  kapustin
  500.  * Code cleanup
  501.  *
  502.  * Revision 1.8  2004/05/17 14:50:56  kapustin
  503.  * Add/remove/rearrange some includes and object declarations
  504.  *
  505.  * Revision 1.7  2004/03/18 16:30:24  grichenk
  506.  * Changed type of seq-align containers from list to vector.
  507.  *
  508.  * Revision 1.6  2003/10/14 19:26:59  kapustin
  509.  * Format void exons properly
  510.  *
  511.  * Revision 1.5  2003/09/26 14:43:18  kapustin
  512.  * Remove exception specifications
  513.  *
  514.  * Revision 1.4  2003/09/15 21:32:03  kapustin
  515.  * Minor code cleanup
  516.  *
  517.  * Revision 1.3  2003/09/12 19:43:04  kapustin
  518.  * Add checking for empty transcript
  519.  *
  520.  * Revision 1.2  2003/09/03 01:19:32  ucko
  521.  * +<iterator> (needed for ostream_iterator<> with some compilers)
  522.  *
  523.  * Revision 1.1  2003/09/02 22:34:49  kapustin
  524.  * Initial revision
  525.  *
  526.  * ===========================================================================
  527.  */