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

生物技术

开发平台:

C/C++

  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: db_blast.cpp,v $
  4.  * PRODUCTION Revision 1000.5  2004/06/01 18:06:04  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.28
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  9. /*  $Id: db_blast.cpp,v 1000.5 2004/06/01 18:06:04 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:  Ilya Dondoshansky
  35.  *
  36.  * File Description:
  37.  *   Database BLAST class interface
  38.  *
  39.  * ===========================================================================
  40.  */
  41. #include <ncbi_pch.hpp>
  42. #include <objmgr/util/sequence.hpp>
  43. #include <objects/seqloc/Seq_loc.hpp>
  44. #include <objects/seqalign/Seq_align_set.hpp>
  45. #include <objects/seqalign/Seq_align.hpp>
  46. #include <algo/blast/api/db_blast.hpp>
  47. #include <algo/blast/api/blast_options.hpp>
  48. #include "blast_seqalign.hpp"
  49. #include "blast_setup.hpp"
  50. // NewBlast includes
  51. #include <algo/blast/core/blast_def.h>
  52. #include <algo/blast/core/blast_util.h>
  53. #include <algo/blast/core/blast_setup.h>
  54. #include <algo/blast/core/lookup_wrap.h>
  55. #include <algo/blast/core/blast_engine.h>
  56. #include <algo/blast/core/blast_message.h>
  57. BEGIN_NCBI_SCOPE
  58. USING_SCOPE(objects);
  59. BEGIN_SCOPE(blast)
  60. void CDbBlast::x_InitFields()
  61. {
  62.     m_ibQuerySetUpDone = false;
  63.     m_ipScoreBlock = NULL;
  64.     m_ipLookupTable = NULL;
  65.     m_ipLookupSegments = NULL;
  66.     m_ipFilteredRegions = NULL;
  67.     m_ipResults = NULL;
  68.     m_ipDiagnostics = NULL;
  69.     m_OptsHandle->SetDbLength(BLASTSeqSrcGetTotLen(m_pSeqSrc));
  70.     m_OptsHandle->SetDbSeqNum(BLASTSeqSrcGetNumSeqs(m_pSeqSrc));
  71. }
  72. CDbBlast::CDbBlast(const TSeqLocVector& queries, BlastSeqSrc* seq_src,
  73.                    EProgram p, RPSInfo* rps_info)
  74.     : m_tQueries(queries), m_pSeqSrc(seq_src), m_pRpsInfo(rps_info) 
  75. {
  76.     m_OptsHandle.Reset(CBlastOptionsFactory::Create(p));
  77.     x_InitFields();
  78. }
  79. CDbBlast::CDbBlast(const TSeqLocVector& queries, BlastSeqSrc* seq_src, 
  80.                    CBlastOptionsHandle& opts, RPSInfo* rps_info)
  81.     : m_tQueries(queries), m_pSeqSrc(seq_src), m_pRpsInfo(rps_info) 
  82. {
  83.     m_OptsHandle.Reset(&opts);    
  84.     x_InitFields();
  85. }
  86. CDbBlast::~CDbBlast()
  87.     x_ResetQueryDs();
  88. }
  89. /// Resets results data structures
  90. void
  91. CDbBlast::x_ResetResultDs()
  92. {
  93.     m_ipResults = Blast_HSPResultsFree(m_ipResults);
  94.     
  95.     m_ipDiagnostics = Blast_DiagnosticsFree(m_ipDiagnostics);
  96.     NON_CONST_ITERATE(TBlastError, itr, m_ivErrors) {
  97.         *itr = Blast_MessageFree(*itr);
  98.     }
  99. }
  100. /// Resets query data structures
  101. void
  102. CDbBlast::x_ResetQueryDs()
  103. {
  104.     m_ibQuerySetUpDone = false;
  105.     // should be changed if derived classes are created
  106.     m_iclsQueries.Reset(NULL);
  107.     m_iclsQueryInfo.Reset(NULL);
  108.     m_ipScoreBlock = BlastScoreBlkFree(m_ipScoreBlock);
  109.     m_ipLookupTable = LookupTableWrapFree(m_ipLookupTable);
  110.     m_ipLookupSegments = ListNodeFreeData(m_ipLookupSegments);
  111.     m_ipFilteredRegions = BlastMaskLocFree(m_ipFilteredRegions);
  112.     x_ResetResultDs();
  113. }
  114. int CDbBlast::SetupSearch()
  115. {
  116.     int status = 0;
  117.     EProgram x_eProgram = m_OptsHandle->GetOptions().GetProgram();
  118.     if ( !m_ibQuerySetUpDone ) {
  119.         double scale_factor;
  120.         x_ResetQueryDs();
  121.         bool translated_query = (x_eProgram == eBlastx || 
  122.                                  x_eProgram == eTblastx);
  123.         SetupQueryInfo(m_tQueries, m_OptsHandle->GetOptions(), 
  124.                        &m_iclsQueryInfo);
  125.         SetupQueries(m_tQueries, m_OptsHandle->GetOptions(), 
  126.                      m_iclsQueryInfo, &m_iclsQueries);
  127.         m_ipScoreBlock = 0;
  128.         if (x_eProgram == eRPSBlast || x_eProgram == eRPSTblastn)
  129.             scale_factor = m_pRpsInfo->aux_info.scale_factor;
  130.         else
  131.             scale_factor = 1.0;
  132.         Blast_Message* blast_message = NULL;
  133.         status = BLAST_MainSetUp(x_eProgram, 
  134.                                  m_OptsHandle->GetOptions().GetQueryOpts(),
  135.                                  m_OptsHandle->GetOptions().GetScoringOpts(),
  136.                                  m_OptsHandle->GetOptions().GetHitSaveOpts(),
  137.                                  m_iclsQueries, m_iclsQueryInfo,
  138.                                  scale_factor,
  139.                                  &m_ipLookupSegments, &m_ipFilteredRegions,
  140.                                  &m_ipScoreBlock, &blast_message);
  141.         if (status != 0) {
  142.             string msg = blast_message ? blast_message->message : 
  143.                 "BLAST_MainSetUp failed";
  144.             Blast_MessageFree(blast_message);
  145.             NCBI_THROW(CBlastException, eInternal, msg.c_str());
  146.         } else if (blast_message) {
  147.             // Non-fatal error message; just save it.
  148.             m_ivErrors.push_back(blast_message);
  149.         }
  150.         
  151.         if (translated_query) {
  152.             /* Filter locations were returned in protein coordinates; 
  153.                convert them back to nucleotide here */
  154.             BlastMaskLocProteinToDNA(&m_ipFilteredRegions, m_tQueries);
  155.         }
  156.         Blast_HSPResultsInit(m_iclsQueryInfo->num_queries, &m_ipResults);
  157.         LookupTableWrapInit(m_iclsQueries, 
  158.             m_OptsHandle->GetOptions().GetLutOpts(), 
  159.             m_ipLookupSegments, m_ipScoreBlock, &m_ipLookupTable, m_pRpsInfo);
  160.         
  161.         m_ibQuerySetUpDone = true;
  162.         m_ipDiagnostics = Blast_DiagnosticsInit();
  163.     }
  164.     return status;
  165. }
  166. void CDbBlast::PartialRun()
  167. {
  168.     SetupSearch();
  169.     m_OptsHandle->GetOptions().Validate();
  170.     RunSearchEngine();
  171. }
  172. TSeqAlignVector
  173. CDbBlast::Run()
  174. {
  175.     m_OptsHandle->GetOptions().Validate();// throws an exception on failure
  176.     SetupSearch();
  177.     RunSearchEngine();
  178.     return x_Results2SeqAlign();
  179. }
  180. /* Comparison function for sorting HSP lists in increasing order of the 
  181.    number of HSPs in a hit. Needed for TrimBlastHSPResults below. */
  182. extern "C" int
  183. compare_hsplist_hspcnt(const void* v1, const void* v2)
  184. {
  185.    BlastHSPList* r1 = *((BlastHSPList**) v1);
  186.    BlastHSPList* r2 = *((BlastHSPList**) v2);
  187.    if (r1->hspcnt < r2->hspcnt)
  188.       return -1;
  189.    else if (r1->hspcnt > r2->hspcnt)
  190.       return 1;
  191.    else
  192.       return 0;
  193. }
  194. /** Removes extra results if a limit is imposed on the total number of HSPs
  195.  * returned. Makes sure that at least 1 HSP is still returned for each
  196.  * database sequence hit. 
  197.  */
  198. void 
  199. CDbBlast::TrimBlastHSPResults()
  200. {
  201.     int total_hsp_limit = m_OptsHandle->GetOptions().GetTotalHspLimit();
  202.     if (total_hsp_limit == 0)
  203.         return;
  204.     Int4 total_hsps = 0;
  205.     Int4 allowed_hsp_num, hsplist_count;
  206.     bool hsp_limit_exceeded = false;
  207.     BlastHitList* hit_list;
  208.     BlastHSPList* hsp_list;
  209.     BlastHSPList** hsplist_array;
  210.     int query_index, subject_index, hsp_index;
  211.     
  212.     for (query_index = 0; query_index < m_ipResults->num_queries; 
  213.          ++query_index) {
  214.         if (!(hit_list = m_ipResults->hitlist_array[query_index]))
  215.             continue;
  216.         /* The count of HSPs is separate for each query. */
  217.         total_hsps = 0;
  218.         hsplist_count = hit_list->hsplist_count;
  219.         hsplist_array = (BlastHSPList**) 
  220.             malloc(hsplist_count*sizeof(BlastHSPList*));
  221.         
  222.         for (subject_index = 0; subject_index < hsplist_count; ++subject_index)
  223.             hsplist_array[subject_index] = 
  224.                 hit_list->hsplist_array[subject_index];
  225.         
  226.         qsort((void*)hsplist_array, hsplist_count,
  227.               sizeof(BlastHSPList*), compare_hsplist_hspcnt);
  228.         
  229.         for (subject_index = 0; subject_index < hsplist_count; 
  230.              ++subject_index) {
  231.             allowed_hsp_num = 
  232.                 ((subject_index+1)*total_hsp_limit)/hsplist_count - total_hsps;
  233.             hsp_list = hsplist_array[subject_index];
  234.             if (hsp_list->hspcnt > allowed_hsp_num) {
  235.                 hsp_limit_exceeded = true;
  236.                 /* Free the extra HSPs */
  237.                 for (hsp_index = allowed_hsp_num; 
  238.                      hsp_index < hsp_list->hspcnt; ++hsp_index)
  239.                     Blast_HSPFree(hsp_list->hsp_array[hsp_index]);
  240.                 hsp_list->hspcnt = allowed_hsp_num;
  241.             }
  242.             total_hsps += hsp_list->hspcnt;
  243.         }
  244.         sfree(hsplist_array);
  245.     }
  246.     if (hsp_limit_exceeded) {
  247.         Blast_Message* blast_message = NULL;
  248.         Blast_MessageWrite(&blast_message, BLAST_SEV_INFO, 0, 0, 
  249.                            "Too many HSPs to save all");
  250.         m_ivErrors.push_back(blast_message);
  251.     }
  252. }
  253. void 
  254. CDbBlast::RunSearchEngine()
  255. {
  256.     if (m_ipLookupTable->lut_type == RPS_LOOKUP_TABLE) {
  257.         BLAST_RPSSearchEngine(m_OptsHandle->GetOptions().GetProgram(), 
  258.             m_iclsQueries, m_iclsQueryInfo, m_pSeqSrc, m_ipScoreBlock,
  259.             m_OptsHandle->GetOptions().GetScoringOpts(), 
  260.             m_ipLookupTable, m_OptsHandle->GetOptions().GetInitWordOpts(), 
  261.             m_OptsHandle->GetOptions().GetExtnOpts(), 
  262.             m_OptsHandle->GetOptions().GetHitSaveOpts(),
  263.             m_OptsHandle->GetOptions().GetEffLenOpts(),
  264.             m_OptsHandle->GetOptions().GetPSIBlastOpts(), 
  265.             m_OptsHandle->GetOptions().GetDbOpts(),
  266.             m_ipResults, m_ipDiagnostics);
  267.     } else {
  268.         BLAST_SearchEngine(m_OptsHandle->GetOptions().GetProgram(),
  269.             m_iclsQueries, m_iclsQueryInfo, m_pSeqSrc, m_ipScoreBlock, 
  270.             m_OptsHandle->GetOptions().GetScoringOpts(), 
  271.             m_ipLookupTable, m_OptsHandle->GetOptions().GetInitWordOpts(), 
  272.             m_OptsHandle->GetOptions().GetExtnOpts(), 
  273.             m_OptsHandle->GetOptions().GetHitSaveOpts(), 
  274.             m_OptsHandle->GetOptions().GetEffLenOpts(), NULL, 
  275.             m_OptsHandle->GetOptions().GetDbOpts(),
  276.             m_ipResults, m_ipDiagnostics);
  277.     }
  278.     m_ipLookupTable = LookupTableWrapFree(m_ipLookupTable);
  279.     m_iclsQueries = BlastSequenceBlkFree(m_iclsQueries);
  280.     /* The following works because the ListNodes' data point to simple
  281.        double-integer structures */
  282.     m_ipLookupSegments = ListNodeFreeData(m_ipLookupSegments);
  283.     /* If a limit is provided for number of HSPs to return, trim the extra
  284.        HSPs here */
  285.     TrimBlastHSPResults();
  286. }
  287. TSeqAlignVector
  288. CDbBlast::x_Results2SeqAlign()
  289. {
  290.     TSeqAlignVector retval;
  291.     retval = BLAST_Results2CSeqAlign(m_ipResults, 
  292.                  m_OptsHandle->GetOptions().GetProgram(),
  293.                  m_tQueries, m_pSeqSrc, 
  294.                  m_OptsHandle->GetOptions().GetScoringOpts(), 
  295.                  m_ipScoreBlock);
  296.     return retval;
  297. }
  298. END_SCOPE(blast)
  299. END_NCBI_SCOPE
  300. /*
  301.  * ===========================================================================
  302.  *
  303.  * $Log: db_blast.cpp,v $
  304.  * Revision 1000.5  2004/06/01 18:06:04  gouriano
  305.  * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.28
  306.  *
  307.  * Revision 1.28  2004/05/21 21:41:02  gorelenk
  308.  * Added PCH ncbi_pch.hpp
  309.  *
  310.  * Revision 1.27  2004/05/17 18:12:29  bealer
  311.  * - Add PSI-Blast support.
  312.  *
  313.  * Revision 1.26  2004/05/14 17:16:12  dondosha
  314.  * BlastReturnStat structure changed to BlastDiagnostics and refactored
  315.  *
  316.  * Revision 1.25  2004/05/07 15:28:41  papadopo
  317.  * add scale factor to BlastMainSetup
  318.  *
  319.  * Revision 1.24  2004/05/05 15:28:56  dondosha
  320.  * Renamed functions in blast_hits.h accordance with new convention Blast_[StructName][Task]
  321.  *
  322.  * Revision 1.23  2004/04/30 16:53:06  dondosha
  323.  * Changed a number of function names to have the same conventional Blast_ prefix
  324.  *
  325.  * Revision 1.22  2004/04/21 17:33:46  madden
  326.  * Rename BlastHSPFree to Blast_HSPFree
  327.  *
  328.  * Revision 1.21  2004/03/24 22:12:46  dondosha
  329.  * Fixed memory leaks
  330.  *
  331.  * Revision 1.20  2004/03/17 14:51:33  camacho
  332.  * Fix compiler errors
  333.  *
  334.  * Revision 1.19  2004/03/16 23:32:28  dondosha
  335.  * Added capability to run RPS BLAST seach; added function x_InitFields; changed mi_ to m_i in member field names
  336.  *
  337.  * Revision 1.18  2004/03/16 14:45:28  dondosha
  338.  * Removed doxygen comments for nonexisting parameters
  339.  *
  340.  * Revision 1.17  2004/03/15 19:57:00  dondosha
  341.  * Merged TwoSequences and Database engines
  342.  *
  343.  * Revision 1.16  2004/03/10 17:37:36  papadopo
  344.  * add (unused) RPS info pointer to LookupTableWrapInit()
  345.  *
  346.  * Revision 1.15  2004/02/24 20:31:39  dondosha
  347.  * Typo fix; removed irrelevant CVS log comments
  348.  *
  349.  * Revision 1.14  2004/02/24 18:19:35  dondosha
  350.  * Removed lookup options argument from call to BLAST_MainSetUp
  351.  *
  352.  * Revision 1.13  2004/02/23 15:45:09  camacho
  353.  * Eliminate compiler warning about qsort
  354.  *
  355.  * Revision 1.12  2004/02/19 21:12:02  dondosha
  356.  * Added handling of error messages; fill info message in TrimBlastHSPResults
  357.  *
  358.  * Revision 1.11  2004/02/18 23:49:08  dondosha
  359.  * Added TrimBlastHSPResults method to remove extra HSPs if limit on total number is provided
  360.  *
  361.  * Revision 1.10  2004/02/13 20:47:20  madden
  362.  * Throw exception rather than ERR_POST if setup fails
  363.  *
  364.  * Revision 1.9  2004/02/13 19:32:55  camacho
  365.  * Removed unnecessary #defines
  366.  *
  367.  * Revision 1.8  2004/01/16 21:51:34  bealer
  368.  * - Changes for Blast4 API
  369.  *
  370.  * Revision 1.7  2003/12/15 23:42:46  dondosha
  371.  * Set database length and number of sequences options in constructor
  372.  *
  373.  * Revision 1.6  2003/12/15 15:56:42  dondosha
  374.  * Added constructor with options handle argument
  375.  *
  376.  * Revision 1.5  2003/12/03 17:41:19  camacho
  377.  * Fix incorrect member initializer list
  378.  *
  379.  * Revision 1.4  2003/12/03 16:45:03  dondosha
  380.  * Use CBlastOptionsHandle class
  381.  *
  382.  * Revision 1.3  2003/11/26 18:36:45  camacho
  383.  * Renaming blast_option*pp -> blast_options*pp
  384.  *
  385.  * Revision 1.2  2003/10/30 21:41:12  dondosha
  386.  * Removed unneeded extra argument from call to BLAST_Results2CSeqAlign
  387.  *
  388.  * Revision 1.1  2003/10/29 22:37:36  dondosha
  389.  * Database BLAST search class methods
  390.  * ===========================================================================
  391.  */