splign_compartment_finder.cpp
上传用户:yhdzpy8989
上传日期:2007-06-13
资源大小:13604k
文件大小:13k
源码类别:

生物技术

开发平台:

C/C++

  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: splign_compartment_finder.cpp,v $
  4.  * PRODUCTION Revision 1000.0  2004/06/01 18:12:33  gouriano
  5.  * PRODUCTION PRODUCTION: IMPORTED [GCC34_MSVC7] Dev-tree R1.5
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  9. /* $Id: splign_compartment_finder.cpp,v 1000.0 2004/06/01 18:12:33 gouriano Exp $
  10. * ===========================================================================
  11. *
  12. *                            PUBLIC DOMAIN NOTICE                          
  13. *               National Center for Biotechnology Information
  14. *                                                                          
  15. *  This software/database is a "United States Government Work" under the   
  16. *  terms of the United States Copyright Act.  It was written as part of    
  17. *  the author's official duties as a United States Government employee and 
  18. *  thus cannot be copyrighted.  This software/database is freely available 
  19. *  to the public for use. The National Library of Medicine and the U.S.    
  20. *  Government have not placed any restriction on its use or reproduction.  
  21. *                                                                          
  22. *  Although all reasonable efforts have been taken to ensure the accuracy  
  23. *  and reliability of the software and data, the NLM and the U.S.          
  24. *  Government do not and cannot warrant the performance or results that    
  25. *  may be obtained by using this software or data. The NLM and the U.S.    
  26. *  Government disclaim all warranties, express or implied, including       
  27. *  warranties of performance, merchantability or fitness for any particular
  28. *  purpose.                                                                
  29. *                                                                          
  30. *  Please cite the author in any work or product based on this material.   
  31. *
  32. * ===========================================================================
  33. *
  34. * Author:  Yuri Kapustin
  35. *
  36. * File Description: Compartment Finder
  37. *                   
  38. * ===========================================================================
  39. */
  40. #include <ncbi_pch.hpp>
  41. #include "splign_compartment_finder.hpp"
  42. #include <corelib/ncbi_limits.hpp>
  43. #include <algo/align/align_exception.hpp>
  44. #include <algorithm>
  45. BEGIN_NCBI_SCOPE
  46. void CCompartmentFinder::CCompartment::UpdateMinMax() {
  47.   m_box[0] = m_box[2] = kMax_UInt;
  48.   m_box[1] = m_box[3] = 0;
  49.   ITERATE(vector<const CHit*>, ph, m_members) {
  50.     
  51.     const CHit* hit = *ph;
  52.     if(size_t(hit->m_ai[0]) < m_box[0]) {
  53.       m_box[0] = hit->m_ai[0];
  54.     }
  55.     if(size_t(hit->m_ai[2]) < m_box[2]) {
  56.       m_box[2] = hit->m_ai[2];
  57.     }
  58.     if(size_t(hit->m_ai[1]) > m_box[1]) {
  59.       m_box[1] = hit->m_ai[1];
  60.     }
  61.     if(size_t(hit->m_ai[3]) > m_box[3]) {
  62.       m_box[3] = hit->m_ai[3];
  63.     }
  64.   }
  65. }
  66. bool CCompartmentFinder::CCompartment::GetStrand() const {
  67.   if(m_members.size()) {
  68.     return m_members[0]->IsStraight();
  69.   }
  70.   else {
  71.         NCBI_THROW( CAlgoAlignException,
  72.                     eInternal,
  73.                     "Strand requested on an empty compartment" );
  74.   }
  75. }
  76. CCompartmentFinder::CCompartmentFinder(vector<CHit>::const_iterator start,
  77.                                        vector<CHit>::const_iterator finish):
  78.   m_intron_min(GetDefaultIntronMin()),
  79.   m_intron_max(GetDefaultIntronMax()),
  80.   m_penalty(GetDefaultPenalty()),
  81.   m_min_coverage(GetDefaultMinCoverage()),
  82.   m_iter(-1)
  83. {
  84.   const size_t dim = finish - start;
  85.   if(dim >  0) {
  86.     m_hits.resize(dim);
  87.     const CHit* p0 = &(*start);
  88.     for(size_t i = 0; i < dim; ++i) {
  89.       m_hits[i] = p0++;
  90.     }
  91.   }
  92. }
  93. size_t CCompartmentFinder::Run()
  94. {
  95.   m_compartments.clear();
  96.   const size_t dimhits = m_hits.size();
  97.   if(dimhits == 0) {
  98.     return 0;
  99.   }
  100.   // sort here to make sure that later at every step
  101.   // all potential sources (hits from where to continue)
  102.   // are visible
  103.   sort(m_hits.begin(), m_hits.end(), CHit::PSubjLow_QueryLow_ptr);
  104.   // insert dummy element
  105.   list<SQueryMark> qmarks; // ordered list of query marks
  106.   qmarks.clear();
  107.   qmarks.push_back(SQueryMark(0, 0, -1));
  108.   const list<SQueryMark>::iterator li_b = qmarks.begin();
  109.   // *** Generic hit iteration (could be optimized of course) ***
  110.   // For every hit:
  111.   // - find out if its query mark already exists (qm_new)
  112.   // - for query mark below the current:
  113.   //   -- check if it could be extended
  114.   //   -- evaluate extension potential
  115.   // - for every other qm:
  116.   //   -- check if its compartment could be terminated
  117.   //   -- evaluate potential of terminating the compartment
  118.   //      to open a new one
  119.   // - select the best potential
  120.   // - if qm_new, create the new query mark, otherwise
  121.   //   update the existing one if the new score is higher
  122.   // - if query mark was created or updated, 
  123.   //   save the hit's status ( previous hit, extension or opening)
  124.   //
  125.   vector<SHitStatus> hitstatus (dimhits);
  126.   for(size_t i = 0; i < dimhits; ++i) {
  127.     const list<SQueryMark>::iterator li_e = qmarks.end();
  128.     const CHit& h = *m_hits[i];
  129.     list<SQueryMark>::iterator li0 = lower_bound(li_b, li_e,
  130.                                                  SQueryMark(h.m_ai[1], 0, -2));
  131.     bool qm_new = (li0 == li_e)? true: (size_t(h.m_ai[1]) < li0->m_coord);
  132.     int best_ext_score = kMin_Int;
  133.     list<SQueryMark>::iterator li_best_ext = li_b;
  134.     int best_open_score = kMin_Int;
  135.     list<SQueryMark>::iterator li_best_open = li_b;
  136.     for( list<SQueryMark>::iterator li = li_b; li != li_e; ++li) {
  137.       const CHit* phc = (li->m_hit == -1)? 0: m_hits[li->m_hit];
  138.       // check if good for extention
  139.       if(li->m_hit < int(i)) {
  140.         size_t q0 = h.m_ai[0], s0 = h.m_ai[2]; // possible continuation point
  141.         bool good;
  142.         if(li->m_hit == -1) {
  143.           good = false;
  144.         }
  145.         else {
  146.           if(phc->m_ai[1] >= h.m_ai[1]) {
  147.             good = false;
  148.           }
  149.           else {
  150.             if(phc->m_ai[1] < h.m_ai[0]) {
  151.               q0 = h.m_ai[0];
  152.               s0 = h.m_ai[2];
  153.             }
  154.             else {
  155.               // find where this point would be on h
  156.               q0 = phc->m_ai[1] + 1;
  157.               s0 += q0 - h.m_ai[0];
  158.             }
  159.             int intron = int(s0) - phc->m_ai[3] - 1;
  160.             good =  int(m_intron_min) <= intron && intron <= int(m_intron_max);
  161.           }          
  162.         }
  163.         if(good) {
  164.           int ext_score = li->m_score +
  165.               int(0.01 * h.m_Idnty * (h.m_ai[1] - q0 + 1));
  166.           if(ext_score > best_ext_score) {
  167.             best_ext_score = ext_score;
  168.             li_best_ext = li;
  169.           }
  170.         }
  171.       }
  172.       
  173.       // check if good for closing and opening
  174.       if(li->m_hit == -1 ||  phc->m_ai[3] < h.m_ai[2]) {
  175.         int score_open = li->m_score - m_penalty +
  176.           int(0.01*h.m_Idnty*(h.m_ai[1] - h.m_ai[0] + 1));
  177.         if(score_open > best_open_score) {
  178.           best_open_score = score_open;
  179.           li_best_open = li;
  180.         }
  181.       }
  182.     }
  183.     SHitStatus::EType hit_type;
  184.     int prev_hit;
  185.     int best_score;
  186.     if(best_ext_score > best_open_score) {
  187.       hit_type = SHitStatus::eExtension;
  188.       prev_hit = li_best_ext->m_hit;
  189.       best_score = best_ext_score;
  190.     }
  191.     else {
  192.       hit_type = SHitStatus::eOpening;
  193.       prev_hit = li_best_open->m_hit;
  194.       best_score = best_open_score;
  195.     }
  196.     bool updated = false;
  197.     if(qm_new) {
  198.       qmarks.insert(li0, SQueryMark(h.m_ai[1], best_score, i));
  199.       updated = true;
  200.     }
  201.     else {
  202.       if(best_score > li0->m_score) {
  203.         li0->m_score = best_score;
  204.         li0->m_hit = i;
  205.         updated = true;
  206.       }
  207.     }
  208.     if(updated) {
  209.       hitstatus[i].m_type = hit_type;
  210.       hitstatus[i].m_prev = prev_hit;
  211.     }
  212.   }
  213.   // *** backtrace ***
  214.   // - find query mark with the highest score
  215.   // - trace it back until the dummy
  216.   list<SQueryMark>::iterator li_best = li_b;
  217.   ++li_best;
  218.   int score_best = kMin_Int;
  219.   for(list<SQueryMark>::iterator li = li_best, li_e = qmarks.end();
  220.       li != li_e; ++li) {
  221.     if(li->m_score > score_best) {
  222.       score_best = li->m_score;
  223.       li_best = li;
  224.     }
  225.   }
  226.   if( int(score_best + m_penalty) >= int(m_min_coverage) ) {
  227.     int i = li_best->m_hit;
  228.     CCompartment* pc = 0;
  229.     bool new_compartment = true;
  230.     while(i != -1) {
  231.       if(new_compartment) {
  232.         new_compartment = false;
  233.         m_compartments.push_back(CCompartment());
  234.         pc = &m_compartments.back();
  235.       }
  236.       pc->AddMember(m_hits[i]);
  237.       
  238.       if(hitstatus[i].m_type == SHitStatus::eOpening) {
  239.         new_compartment = true;
  240.       }
  241.       i = hitstatus[i].m_prev;
  242.     }
  243.   }
  244.   return m_compartments.size();
  245. }
  246. CCompartmentFinder::CCompartment* CCompartmentFinder::GetFirst()
  247. {
  248.   if(m_compartments.size()) {
  249.     m_iter = 0;
  250.     return &m_compartments[m_iter++];
  251.   }
  252.   else {
  253.     m_iter = -1;
  254.     return 0;
  255.   }
  256. }
  257. void CCompartmentFinder::OrderCompartments(void) {
  258.   
  259.   for(int i = 0, dim = m_compartments.size(); i < dim; ++i) {
  260.     m_compartments[i].UpdateMinMax();
  261.   }
  262.   sort(m_compartments.begin(), m_compartments.end(), PLowerSubj);
  263. }
  264. CCompartmentFinder::CCompartment* CCompartmentFinder::GetNext()
  265. {
  266.   const size_t dim = m_compartments.size();
  267.   if(m_iter == -1 || m_iter >= int(dim)) {
  268.     m_iter = -1;
  269.     return 0;
  270.   }
  271.   else {
  272.     return &m_compartments[m_iter++];
  273.   }
  274. }
  275. CCompartmentAccessor::CCompartmentAccessor(vector<CHit>::iterator istart,
  276.                                            vector<CHit>::iterator ifinish,
  277.                                            size_t comp_penalty,
  278.                                            size_t min_coverage)
  279. {
  280.     // separate strands for CompartmentFinder
  281.     vector<CHit>::iterator ib = istart, ie = ifinish, ii = ib, iplus_beg = ie;
  282.     sort(ib, ie, CHit::PPrecedeStrand);
  283.     size_t minus_subj_min = kMax_UInt, minus_subj_max = 0;
  284.     for(ii = ib; ii != ie; ++ii) {
  285.       if(ii->IsStraight()) {
  286.         iplus_beg = ii;
  287.         break;
  288.       }
  289.       else {
  290.         if(size_t(ii->m_ai[2]) < minus_subj_min) {
  291.           minus_subj_min = ii->m_ai[2];
  292.         }
  293.         if(size_t(ii->m_ai[3]) > minus_subj_max) {
  294.           minus_subj_max = ii->m_ai[3];
  295.         }
  296.       }
  297.     }
  298.     // minus
  299.     {{
  300.       // flip
  301.       for(ii = ib; ii != iplus_beg; ++ii) {
  302.         int s0 = minus_subj_min + minus_subj_max - ii->m_ai[3];
  303.         int s1 = minus_subj_min + minus_subj_max - ii->m_ai[2];
  304.         ii->m_an[2] = ii->m_ai[2] = s0;
  305.         ii->m_an[3] = ii->m_ai[3] = s1;
  306.       }
  307.       CCompartmentFinder finder (ib, iplus_beg);
  308.       finder.SetMinCoverage(min_coverage);
  309.       finder.SetPenalty(comp_penalty);
  310.       finder.Run();
  311.       // un-flip
  312.       for(ii = ib; ii != iplus_beg; ++ii) {
  313.         int s0 = minus_subj_min + minus_subj_max - ii->m_ai[3];
  314.         int s1 = minus_subj_min + minus_subj_max - ii->m_ai[2];
  315.         ii->m_an[3] = ii->m_ai[2] = s0;
  316.         ii->m_an[2] = ii->m_ai[3] = s1;
  317.       }
  318.       x_Copy2Pending(finder);
  319.     }}
  320.     // plus
  321.     {{
  322.       CCompartmentFinder finder (iplus_beg, ie);
  323.       finder.SetMinCoverage(min_coverage);
  324.       finder.SetPenalty(comp_penalty);
  325.       finder.Run();
  326.       x_Copy2Pending(finder);
  327.     }}
  328. }
  329. void CCompartmentAccessor::x_Copy2Pending(CCompartmentFinder& finder)
  330. {
  331.   finder.OrderCompartments();
  332.   // copy to pending
  333.   for(CCompartmentFinder::CCompartment* compartment = finder.GetFirst(); compartment;
  334.       compartment = finder.GetNext()) {
  335.     
  336.     m_pending.push_back(vector<CHit>(0));
  337.     vector<CHit>& vh = m_pending.back();
  338.     
  339.     for(const CHit* ph = compartment->GetFirst(); ph; ph = compartment->GetNext()) {
  340.       vh.push_back(*ph);
  341.     }
  342.     
  343.     const size_t* box = compartment->GetBox();
  344.     m_ranges.push_back(box[0]-1);
  345.     m_ranges.push_back(box[1]-1);
  346.     m_ranges.push_back(box[2]-1);
  347.     m_ranges.push_back(box[3]-1);
  348.     m_strands.push_back(compartment->GetStrand());
  349.   }
  350. }
  351. bool CCompartmentAccessor::GetFirst(vector<CHit>& compartment) {
  352.   m_iter = 0;
  353.   return GetNext(compartment);
  354. }
  355. bool CCompartmentAccessor::GetNext(vector<CHit>& compartment) {
  356.   compartment.clear();
  357.   if(m_iter < m_pending.size()) {
  358.     compartment = m_pending[m_iter++];
  359.     return true;
  360.   }
  361.   else {
  362.     return false;
  363.   }
  364. }
  365. END_NCBI_SCOPE
  366. /*
  367.  * ===========================================================================
  368.  * $Log: splign_compartment_finder.cpp,v $
  369.  * Revision 1000.0  2004/06/01 18:12:33  gouriano
  370.  * PRODUCTION: IMPORTED [GCC34_MSVC7] Dev-tree R1.5
  371.  *
  372.  * Revision 1.5  2004/05/24 16:13:57  gorelenk
  373.  * Added PCH ncbi_pch.hpp
  374.  *
  375.  * Revision 1.4  2004/05/18 21:43:40  kapustin
  376.  * Code cleanup
  377.  *
  378.  * Revision 1.3  2004/04/23 14:37:44  kapustin
  379.  * *** empty log message ***
  380.  *
  381.  * ==========================================================================
  382.  */