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

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

生物技术

开发平台：

C/C++

blast_app.cpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: blast_app.cpp,v $
  4.  * PRODUCTION Revision 1000.8  2004/06/01 18:06:35  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.49
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. static char const rcsid[] = "$Id: blast_app.cpp,v 1000.8 2004/06/01 18:06:35 gouriano Exp $";
  12. /*
  13. * ===========================================================================
  14. *
  15. *                            PUBLIC DOMAIN NOTICE
  16. *               National Center for Biotechnology Information
  17. *
  18. *  This software/database is a "United States Government Work" under the
  19. *  terms of the United States Copyright Act.  It was written as part of
  20. *  the author's offical duties as a United States Government employee and
  21. *  thus cannot be copyrighted.  This software/database is freely available
  22. *  to the public for use. The National Library of Medicine and the U.S.
  23. *  Government have not placed any restriction on its use or reproduction.
  24. *
  25. *  Although all reasonable efforts have been taken to ensure the accuracy
  26. *  and reliability of the software and data, the NLM and the U.S.
  27. *  Government do not and cannot warrant the performance or results that
  28. *  may be obtained by using this software or data. The NLM and the U.S.
  29. *  Government disclaim all warranties, express or implied, including
  30. *  warranties of performance, merchantability or fitness for any particular
  31. *  purpose.
  32. *
  33. *  Please cite the author in any work or product based on this material.
  34. *
  35. * ===========================================================================*/
  37. /*****************************************************************************
  39. File name: blast_app.cpp
  41. Author: Ilya Dondoshansky
  43. Contents: C++ driver for running BLAST
  45. ******************************************************************************/
  47. #include <ncbi_pch.hpp>
  48. #include <corelib/ncbiapp.hpp>
  49. #include <corelib/ncbifile.hpp>
  51. #include <objects/seq/seq__.hpp>
  52. #include <objects/seq/seqport_util.hpp>
  53. #include <objects/seqfeat/Genetic_code_table.hpp>
  54. #include <objmgr/object_manager.hpp>
  55. #include <objtools/data_loaders/blastdb/bdbloader.hpp>
  56. #include <objtools/data_loaders/genbank/gbloader.hpp>
  57. #include <objmgr/util/sequence.hpp>
  59. #include <algo/blast/api/blast_options.hpp>
  60. #include <algo/blast/api/blast_nucl_options.hpp>
  61. #include <algo/blast/api/db_blast.hpp>
  62. #include <algo/blast/api/blast_aux.hpp>
  63. #include "blast_input.hpp"
  64. #include <algo/blast/api/seqdb_src.hpp>
  65. #include <objtools/alnmgr/util/blast_format.hpp>
  67. #ifndef NCBI_C_TOOLKIT
  68. #define NCBI_C_TOOLKIT
  69. #endif
  71. // C include files
  72. #include <algo/blast/core/blast_setup.h>
  73. #include <algo/blast/core/blast_util.h>
  74. #include <algo/blast/core/lookup_wrap.h>
  75. #include <algo/blast/core/blast_engine.h>
  77. // For writing out seqalign only
  78. #include <ctools/asn_converter.hpp>
  79. #include <objalign.h>
  81. USING_NCBI_SCOPE;
  82. USING_SCOPE(blast);
  83. USING_SCOPE(objects);
  85. class CBlastApplication : public CNcbiApplication
  86. {
  87. private:
  88.     virtual void Init(void);
  89.     virtual int Run(void);
  90.     virtual void Exit(void);
  91.     void InitScope(void);
  92.     void InitOptions(void);
  93.     void SetOptions(const CArgs& args);
  94.     void ProcessCommandLineArgs(CBlastOptionsHandle* opt, 
  95.                                 BlastSeqSrc* bssp, RPSInfo *rps_info);
  96.     int BlastSearch(void);
  97.     void RegisterBlastDbLoader(char* dbname, bool is_na);
  98.     void FormatResults(CDbBlast& blaster, TSeqAlignVector& seqalignv);
  99.     CRef<CObjectManager> m_ObjMgr;
  100.     CRef<CScope>         m_Scope;
  101. };
  103. void CBlastApplication::Init(void)
  104. {
  105.     HideStdArgs(fHideLogfile | fHideConffile | fHideVersion);
  106.     auto_ptr<CArgDescriptions> arg_desc(new CArgDescriptions);
  108.     arg_desc->SetUsageContext(GetArguments().GetProgramBasename(), "basic local alignment search tool");
  110.     arg_desc->AddKey("program", "program", "Type of BLAST program",
  111.         CArgDescriptions::eString);
  112.     arg_desc->SetConstraint
  113.         ("program", &(*new CArgAllow_Strings, 
  114.                 "blastp", "blastn", "blastx", "tblastn", "tblastx", 
  115.                 "rpsblast", "rpstblastn"));
  116.     arg_desc->AddDefaultKey("query", "query", "Query file name",
  117.                      CArgDescriptions::eInputFile, "stdin");
  118.     arg_desc->AddKey("db", "database", "BLAST database name",
  119.                      CArgDescriptions::eString);
  120.     arg_desc->AddDefaultKey("strand", "strand", 
  121.         "Query strands to search: 1 forward, 2 reverse, 0,3 both",
  122.         CArgDescriptions::eInteger, "0");
  123.     arg_desc->SetConstraint("strand", new CArgAllow_Integers(0,3));
  125.     arg_desc->AddDefaultKey("filter", "filter", "Filtering option",
  126.                             CArgDescriptions::eString, "T");
  127.     arg_desc->AddDefaultKey("lcase", "lcase", "Should lower case be masked?",
  128.                             CArgDescriptions::eBoolean, "F");
  129.     arg_desc->AddDefaultKey("lookup", "lookup", 
  130.         "Type of lookup table: 0 default, 1 megablast",
  131.         CArgDescriptions::eInteger, "0");
  132.     arg_desc->AddDefaultKey("matrix", "matrix", "Scoring matrix name",
  133.                             CArgDescriptions::eString, "BLOSUM62");
  134.     arg_desc->AddDefaultKey("mismatch", "penalty", "Penalty score for a mismatch",
  135.                             CArgDescriptions::eInteger, "0");
  136.     arg_desc->AddDefaultKey("match", "reward", "Reward score for a match",
  137.                             CArgDescriptions::eInteger, "0");
  138.     arg_desc->AddDefaultKey("word", "wordsize", "Word size",
  139.                             CArgDescriptions::eInteger, "0");
  140.     arg_desc->AddDefaultKey("templen", "templen", 
  141.         "Discontiguous word template length",
  142.         CArgDescriptions::eInteger, "0");
  143.     arg_desc->SetConstraint("templen", 
  144.                             &(*new CArgAllow_Strings, "0", "16", "18", "21"));
  146.     arg_desc->AddDefaultKey("templtype", "templtype", 
  147.         "Discontiguous word template type",
  148.         CArgDescriptions::eInteger, "0");
  149.     arg_desc->SetConstraint("templtype", new CArgAllow_Integers(0,2));
  151.     arg_desc->AddDefaultKey("thresh", "threshold", 
  152.         "Score threshold for neighboring words",
  153.         CArgDescriptions::eInteger, "0");
  154.     arg_desc->AddDefaultKey("window","window", "Window size for two-hit extension",
  155.                             CArgDescriptions::eInteger, "0");
  156.     arg_desc->AddDefaultKey("scantype", "scantype", 
  157.         "Method for scanning the database: 0 traditional, 1 AG",
  158.         CArgDescriptions::eInteger, "1");
  159.     arg_desc->SetConstraint("scantype", new CArgAllow_Integers(0,1));
  161.     arg_desc->AddDefaultKey("varword", "varword", 
  162.         "Should variable word size be used?",
  163.         CArgDescriptions::eBoolean, "F");
  164.     arg_desc->AddDefaultKey("stride","stride", "Database scanning stride",
  165.                             CArgDescriptions::eInteger, "0");
  166.     arg_desc->AddDefaultKey("xungap", "xungapped", 
  167.         "X-dropoff value for ungapped extensions",
  168.         CArgDescriptions::eDouble, "0");
  169.     arg_desc->AddDefaultKey("ungapped", "ungapped", 
  170.         "Perform only an ungapped alignment search?",
  171.         CArgDescriptions::eBoolean, "F");
  172.     arg_desc->AddDefaultKey("greedy", "greedy", 
  173.         "Use greedy algorithm for gapped extensions: 0 no, 1 one-step, 2 two-step, 3 two-step with ungapped",
  174.         CArgDescriptions::eInteger, "0");
  175.     arg_desc->AddDefaultKey("gopen", "gapopen", "Penalty for opening a gap",
  176.                             CArgDescriptions::eInteger, "0");
  177.     arg_desc->AddDefaultKey("gext", "gapext", "Penalty for extending a gap",
  178.                             CArgDescriptions::eInteger, "0");
  179.     arg_desc->AddDefaultKey("xgap", "xdrop", 
  180.         "X-dropoff value for preliminary gapped extensions",
  181.         CArgDescriptions::eDouble, "0");
  182.     arg_desc->AddDefaultKey("xfinal", "xfinal", 
  183.         "X-dropoff value for final gapped extensions with traceback",
  184.         CArgDescriptions::eDouble, "0");
  185.     arg_desc->AddDefaultKey("evalue", "evalue", 
  186.         "E-value threshold for saving hits",
  187.         CArgDescriptions::eDouble, "0");
  188.     arg_desc->AddDefaultKey("searchsp", "searchsp", 
  189.         "Virtual search space to be used for statistical calculations",
  190.         CArgDescriptions::eDouble, "0");
  191.     arg_desc->AddDefaultKey("perc", "percident", 
  192.         "Percentage of identities cutoff for saving hits",
  193.         CArgDescriptions::eDouble, "0");
  194.     arg_desc->AddDefaultKey("hitlist", "hitlist_size",
  195.         "How many best matching sequences to find?",
  196.         CArgDescriptions::eInteger, "500");
  197.     arg_desc->AddDefaultKey("descr", "descriptions",
  198.         "How many matching sequence descriptions to show?",
  199.         CArgDescriptions::eInteger, "500");
  200.     arg_desc->AddDefaultKey("align", "alignments", 
  201.         "How many matching sequence alignments to show?",
  202.         CArgDescriptions::eInteger, "250");
  203.     arg_desc->AddOptionalKey("out", "outfile", 
  204.         "File name for writing output",
  205.         CArgDescriptions::eOutputFile);
  206.     arg_desc->AddDefaultKey("format", "format", 
  207.         "How to format the results?",
  208.         CArgDescriptions::eInteger, "0");
  209.     arg_desc->AddDefaultKey("html", "html", "Produce HTML output?",
  210.                             CArgDescriptions::eBoolean, "F");
  211.     arg_desc->AddDefaultKey("gencode", "gencode", "Query genetic code",
  212.                             CArgDescriptions::eInteger, "1");
  213.     arg_desc->AddDefaultKey("dbgencode", "dbgencode", "Database genetic code",
  214.                             CArgDescriptions::eInteger, "1");
  215.     arg_desc->AddDefaultKey("maxintron", "maxintron", 
  216.                             "Longest allowed intron length for linking HSPs",
  217.                             CArgDescriptions::eInteger, "0");
  218.     arg_desc->AddDefaultKey("frameshift", "frameshift",
  219.                             "Frame shift penalty (blastx only)",
  220.                             CArgDescriptions::eInteger, "0");
  221.     arg_desc->AddOptionalKey("asnout", "seqalignasn", 
  222.         "File name for writing the seqalign results in ASN.1 form",
  223.         CArgDescriptions::eOutputFile);
  225.     arg_desc->AddDefaultKey("qstart", "query_start",
  226.                             "Starting offset in query location",
  227.                             CArgDescriptions::eInteger, "0");
  228.     arg_desc->AddDefaultKey("qend", "query_end",
  229.                             "Ending offset in query location",
  230.                             CArgDescriptions::eInteger, "0");
  232.     arg_desc->AddOptionalKey("pattern", "phipattern",
  233.                             "Pattern for PHI-BLAST",
  234.                              CArgDescriptions::eString);
  236.     arg_desc->AddOptionalKey("pssm", "pssm", 
  237.         "File name for uploading a PSSM for PSI BLAST seach",
  238.         CArgDescriptions::eInputFile);
  240.     arg_desc->AddOptionalKey("dbrange", "databaserange",
  241.                             "Range of ordinal ids in the BLAST database.n"
  242.                              "Format: "oid1 oid2"",
  243.                              CArgDescriptions::eString);
  245.     SetupArgDescriptions(arg_desc.release());
  246. }
  248. void 
  249. CBlastApplication::InitScope(void)
  250. {
  251.     if (m_Scope.Empty()) {
  252.         m_ObjMgr.Reset(new CObjectManager());
  253.         m_ObjMgr->RegisterDataLoader(*new CGBDataLoader("ID", 0, 2),
  254.                 CObjectManager::eDefault);
  256.         m_Scope.Reset(new CScope(*m_ObjMgr));
  257.         m_Scope->AddDefaults();
  258.         _TRACE("Blast2seqApp: Initializing scope");
  259.     }
  260. }
  262. void 
  263. CBlastApplication::RegisterBlastDbLoader(char *dbname, bool db_is_na)
  264. {
  265.     m_ObjMgr->RegisterDataLoader(*new CBlastDbDataLoader("BLASTDB", dbname, 
  266.                   db_is_na? (CBlastDbDataLoader::eNucleotide) : 
  267.                             (CBlastDbDataLoader::eProtein)),  
  268.                   CObjectManager::eDefault);
  269. }
  271. void
  272. CBlastApplication::ProcessCommandLineArgs(CBlastOptionsHandle* opts_handle, 
  273.                                           BlastSeqSrc* bssp, RPSInfo *rps_info)
  274. {
  275.     CArgs args = GetArgs();
  276.     CBlastOptions& opt = opts_handle->SetOptions();
  277.     EProgram program_number = opt.GetProgram();
  279.     if (args["strand"].AsInteger()) {
  280.         switch (args["strand"].AsInteger()) {
  281.         case 1: opt.SetStrandOption(eNa_strand_plus); break;
  282.         case 2: opt.SetStrandOption(eNa_strand_minus); break;
  283.         case 3: opt.SetStrandOption(eNa_strand_both); break;
  284.         default: abort();
  285.         }
  286.     }
  288.     opt.SetFilterString(args["filter"].AsString().c_str());
  289.     // FIXME: Handle lcase masking
  291.     // If lookup table type argument value is 0, the type will be set correctly
  292.     // automatically. Value 1 corresponds to megablast lookup table;
  293.     switch (args["lookup"].AsInteger()) {
  294.     case 1:
  295.         opt.SetLookupTableType(MB_LOOKUP_TABLE);
  296.         break;
  297.     default:
  298.         break;
  299.     }
  301.     if (program_number == eRPSBlast ||
  302.         program_number == eRPSTblastn) {
  303.         ASSERT(rps_info != NULL);
  304.         opt.SetGapOpeningCost(rps_info->aux_info.gap_open_penalty);
  305.         opt.SetGapExtensionCost(rps_info->aux_info.gap_extend_penalty);
  306.     }
  307.     else {
  308.         if (args["matrix"]) {
  309.             opt.SetMatrixName(args["matrix"].AsString().c_str());
  310.         }
  311.         if (args["gopen"].AsInteger() || args["greedy"].AsInteger()) {
  312.             opt.SetGapOpeningCost(args["gopen"].AsInteger());
  313.         }
  314.         if (args["gext"].AsInteger() || args["greedy"].AsInteger()) {
  315.             opt.SetGapExtensionCost(args["gext"].AsInteger());
  316.         }
  317.     }
  319.     if (args["mismatch"].AsInteger()) {
  320.         opt.SetMismatchPenalty(args["mismatch"].AsInteger());
  321.     }
  322.     if (args["match"].AsInteger()) {
  323.         opt.SetMatchReward(args["match"].AsInteger());
  324.     }
  325.     if (args["word"].AsInteger()) {
  326.         opt.SetWordSize(args["word"].AsInteger());
  327.     }
  328.     if (args["templen"].AsInteger()) {
  329.         opt.SetMBTemplateLength(args["templen"].AsInteger());
  330.     }
  331.     if (args["templtype"].AsInteger()) {
  332.         opt.SetMBTemplateType(args["templtype"].AsInteger());
  333.     }
  334.     if (args["thresh"].AsInteger()) {
  335.         opt.SetWordThreshold(args["thresh"].AsInteger());
  336.     }
  337.     if (args["window"].AsInteger()) {
  338.         opt.SetWindowSize(args["window"].AsInteger());
  339.     }
  341.     // The next 3 apply to nucleotide searches only
  342.     string program = args["program"].AsString();
  343.     if (program == "blastn") {
  344.         // Setting seed extension method involves changing the scanning 
  345.         // stride as well, which is handled in the derived 
  346.         // CBlastNucleotideOptionsHandle class, but not in the base
  347.         // CBlastOptionsHandle class.
  348.         CBlastNucleotideOptionsHandle* nucl_handle = 
  349.             dynamic_cast<CBlastNucleotideOptionsHandle*>(opts_handle);
  350.         if (!args["templen"].AsInteger()) {
  351.             opt.SetVariableWordSize(args["varword"].AsBoolean());
  352.             switch(args["scantype"].AsInteger()) {
  353.             case 1:
  354.                 nucl_handle->SetSeedExtensionMethod(eRightAndLeft);
  355.                 break;
  356.             default:
  357.                 nucl_handle->SetSeedExtensionMethod(eRight);
  358.                 break;
  359.             }
  360.         } else {
  361.             // Discontiguous Mega BLAST: only one extension method.
  362.             nucl_handle->SetSeedExtensionMethod(eRight); 
  363.         }
  364.         // Override the scan step value if it is set by user
  365.         if (args["stride"].AsInteger()) {
  366.             opt.SetScanStep(args["stride"].AsInteger());
  367.         }
  368.     }
  370.     if (args["xungap"].AsDouble()) {
  371.         opt.SetXDropoff(args["xungap"].AsDouble());
  372.     }
  374.     if (args["ungapped"].AsBoolean()) {
  375.         opt.SetGappedMode(false);
  376.     }
  378.     switch (args["greedy"].AsInteger()) {
  379.     case 1: /* Immediate greedy gapped extension with traceback */
  380.         opt.SetGapExtnAlgorithm(eGreedyWithTracebackExt);
  381.         opt.SetUngappedExtension(false);
  382.         break;
  383.     case 2: /* Two-step greedy extension, no ungapped extension */
  384.         opt.SetGapExtnAlgorithm(eGreedyExt);
  385.         opt.SetUngappedExtension(false);
  386.         break;
  387.     case 3: /* Two-step greedy extension after ungapped extension*/
  388.         opt.SetGapExtnAlgorithm(eGreedyExt);
  389.         opt.SetUngappedExtension(true);
  390.         break;
  391.     default: break;
  392.     }
  395.     if (args["xgap"].AsDouble()) {
  396.         opt.SetGapXDropoff(args["xgap"].AsDouble());
  397.     }
  398.     if (args["xfinal"].AsDouble()) {
  399.         opt.SetGapXDropoffFinal(args["xfinal"].AsDouble());
  400.     }
  402.     if (args["evalue"].AsDouble()) {
  403.         opt.SetEvalueThreshold(args["evalue"].AsDouble());
  404.     }
  406.     if (args["searchsp"].AsDouble()) {
  407.         opt.SetEffectiveSearchSpace((Int8) args["searchsp"].AsDouble());
  408.     } else if (bssp) {
  409.         opt.SetDbLength(BLASTSeqSrcGetTotLen(bssp));
  410.         opt.SetDbSeqNum(BLASTSeqSrcGetNumSeqs(bssp));
  411.     }
  413.     if (args["hitlist"].AsInteger()) {
  414.         opt.SetHitlistSize(args["hitlist"].AsInteger());
  415.         /* Hitlist size for preliminary alignments is increased, unless 
  416.            no traceback is performed. */
  417.         if (args["ungapped"].AsBoolean() || args["greedy"].AsInteger() == 1) {
  418.             opt.SetPrelimHitlistSize(args["hitlist"].AsInteger());
  419.         } else {
  420.             opt.SetPrelimHitlistSize(MIN(2*args["hitlist"].AsInteger(),
  421.                                          args["hitlist"].AsInteger() + 50));
  422.         }
  423.     }
  425.     if (args["perc"].AsDouble()) {
  426.         opt.SetPercentIdentity(args["perc"].AsDouble());
  427.     }
  429.     if (args["gencode"].AsInteger()) {
  430.         opt.SetQueryGeneticCode(args["gencode"].AsInteger());
  431.     }
  432.   
  433.     if ((program == "tblastn" || program == "tblastx") &&
  434.         args["dbgencode"].AsInteger() != BLAST_GENETIC_CODE) {
  435.         opt.SetDbGeneticCode(args["dbgencode"].AsInteger());
  436.     }
  438.     if (args["maxintron"].AsInteger()) {
  439.         opt.SetLongestIntronLength(args["maxintron"].AsInteger());
  440.     }
  441.     if (args["frameshift"].AsInteger()) {
  442.         opt.SetFrameShiftPenalty(args["frameshift"].AsInteger());
  443.         opt.SetOutOfFrameMode();
  444.     }
  446.     if (args["pattern"]) {
  447.         opt.SetPHIPattern(args["pattern"].AsString().c_str(),
  448.                                  (program == "blastn"));
  449.     }
  451.     if (args["pssm"]) {
  452. #if 0
  453.         FILE* pssmfp = fopen(args["pssm"].AsInputFile(), "r");
  454.         // Read the pssm string and fill the posMatrix
  455.         bool use_pssm = TRUE;
  456.         fclose(pssmfp);
  457. #endif
  458.     }
  460.     return;
  461. }
  463. static CDisplaySeqalign::TranslatedFrames 
  464. Context2TranslatedFrame(int context)
  465. {
  466.     switch (context) {
  467.     case 1: return CDisplaySeqalign::ePlusStrand1;
  468.     case 2: return CDisplaySeqalign::ePlusStrand2;
  469.     case 3: return CDisplaySeqalign::ePlusStrand3;
  470.     case 4: return CDisplaySeqalign::eMinusStrand1;
  471.     case 5: return CDisplaySeqalign::eMinusStrand2;
  472.     case 6: return CDisplaySeqalign::eMinusStrand3;
  473.     default: return CDisplaySeqalign::eFrameNotSet;
  474.     }
  475. }
  477. #define NUM_FRAMES 6
  479. static TSeqLocInfoVector
  480. BlastMaskLoc2CSeqLoc(const BlastMaskLoc* mask, const TSeqLocVector& slp,
  481.     EProgram program)
  482. {
  483.     TSeqLocInfoVector retval;
  484.     int frame, num_frames;
  485.     bool translated_query;
  486.     int index;
  488.     translated_query = (program == eBlastx ||
  489.                         program == eTblastx);
  491.     num_frames = (translated_query ? NUM_FRAMES : 1);
  493.     TSeqLocInfo mask_info_list;
  495.     for (index = 0; index < (int)slp.size(); ++index) {
  496.         mask_info_list.clear();
  498.         if (!mask) {
  499.             retval.push_back(mask_info_list);
  500.             continue;
  501.         }
  502.         for ( ; mask && mask->index < index*num_frames;
  503.               mask = mask->next);
  504.         BlastSeqLoc* loc;
  505.         CDisplaySeqalign::SeqlocInfo* seqloc_info =
  506.             new CDisplaySeqalign::SeqlocInfo;
  508.         for ( ; mask && mask->index < (index+1)*num_frames;
  509.               mask = mask->next) {
  510.             frame = (translated_query ? (mask->index % num_frames) + 1 : 0);
  513.             for (loc = mask->loc_list; loc; loc = loc->next) {
  514.                 seqloc_info->frame = Context2TranslatedFrame(frame);
  515.                 CRef<CSeq_loc> seqloc(new CSeq_loc());
  516.                 seqloc->SetInt().SetFrom(((SSeqRange*) loc->ptr)->left);
  517.                 seqloc->SetInt().SetTo(((SSeqRange*) loc->ptr)->right);
  518.                 seqloc->SetInt().SetId(*(const_cast<CSeq_id*>(&sequence::GetId(*
  519. slp[index].seqloc, slp[index].scope))));
  521.                 seqloc_info->seqloc = seqloc;
  522.                 mask_info_list.push_back(seqloc_info);
  523.             }
  524.         }
  525.         retval.push_back(mask_info_list);
  526.     }
  528.     return retval;
  529. }
  531. void CBlastApplication::FormatResults(CDbBlast& blaster, 
  532.                                       TSeqAlignVector& seqalignv)
  533. {
  534.     CArgs args = GetArgs();
  536.     if (args["asnout"]) {
  537.         auto_ptr<CObjectOStream> asnout(
  538.             CObjectOStream::Open(args["asnout"].AsString(), eSerial_AsnText));
  539.         unsigned int query_index;
  540.         for (query_index = 0; query_index < seqalignv.size(); ++query_index)
  541.         {
  542.             if (!seqalignv[query_index]->IsSet())
  543.                 continue;
  544.             *asnout << *seqalignv[query_index];
  545.         }
  546.     }
  548.     if (args["out"]) {
  549.         EProgram program = blaster.SetOptions().GetProgram();
  551.         char* dbname = const_cast<char*>(args["db"].AsString().c_str());
  552.         bool db_is_na = (program == eBlastn || program == eTblastn || 
  553.                          program == eTblastx);
  555.         /* Revert RPS program names to their conventional
  556.            counterparts, to avoid confusing the C toolkit
  557.            formatter */
  559.         if (program == eRPSBlast)
  560.             program = eBlastp;
  561.         if (program == eRPSTblastn)
  562.             program = eTblastn;
  564.         CBlastFormatOptions* format_options = 
  565.             new CBlastFormatOptions(program, args["out"].AsOutputFile());
  566.         
  567.         format_options->SetAlignments(args["align"].AsInteger());
  568.         format_options->SetDescriptions(args["descr"].AsInteger());
  569.         format_options->SetAlignView(args["format"].AsInteger());
  570.         format_options->SetHtml(args["html"].AsBoolean());
  571.         
  572. #ifdef C_FORMATTING
  573.         if (dbname) {
  574.             BLAST_PrintOutputHeader(format_options, 
  575.                 args["greedy"].AsBoolean(), dbname, db_is_na);
  576.         }
  577. #endif
  578.         
  579.         RegisterBlastDbLoader(dbname, db_is_na);
  580.         /* Format the results */
  581.         TSeqLocInfoVector maskv =
  582.             BlastMaskLoc2CSeqLoc(blaster.GetFilteredQueryRegions(), 
  583.                               blaster.GetQueries(), program);
  584.         
  585.         if (BLAST_FormatResults(seqalignv, program, blaster.GetQueries(), 
  586.                 maskv, format_options, blaster.GetOptions().GetOutOfFrameMode())) {
  587.             NCBI_THROW(CBlastException, eInternal, 
  588.                        "Error in formatting results");
  589.         }
  590.         
  591. #ifdef C_FORMATTING
  592.         PrintOutputFooter(program, format_options, score_options, 
  593.             m_sbp, lookup_options, word_options, ext_options, hit_options, 
  594.             blaster.GetQueryInfo(), dbname, blaster.GetDiagnostics());
  595. #endif
  596.         
  597.     }
  598. }
  600. static Int2 BLAST_FillRPSInfo( RPSInfo **ppinfo, Nlm_MemMap **rps_mmap,
  601.                                Nlm_MemMap **rps_pssm_mmap, const char* dbname )
  602. {
  603.    char pathname[PATH_MAX];
  604.    char filename[PATH_MAX];
  605.    RPSInfo *info;
  606.    FILE *auxfile;
  607.    Int4 i;
  608.    Int4 seq_size;
  609.    Int4 num_db_seqs;
  610.    Nlm_MemMapPtr lut_mmap;
  611.    Nlm_MemMapPtr pssm_mmap;
  612.    char buffer[PATH_MAX];
  614.    info = (RPSInfo *)malloc(sizeof(RPSInfo));
  615.    if (info == NULL)
  616.       ErrPostEx(SEV_FATAL, 1, 0, "Memory allocation failed");
  618.    /* construct the full path to the DB file. Look in
  619.       the local directory, then BLASTDB environment 
  620.       variable (if any), then .ncbirc */
  622.    sprintf(filename, "%s.loo", dbname);
  624.    if (FileLength(filename) > 0) {
  625.       strcpy(pathname, dbname);
  626.    } else {
  627. #ifdef OS_UNIX
  628.       if (getenv("BLASTDB"))
  629.          Nlm_GetAppParam("NCBI", "BLAST", "BLASTDB", 
  630.                          getenv("BLASTDB"), pathname, PATH_MAX);
  631.       else
  632. #endif
  633.          Nlm_GetAppParam ("NCBI", "BLAST", "BLASTDB", 
  634.                           BLASTDB_DIR, pathname, PATH_MAX);
  635.       sprintf(filename, "%s%s%s", pathname, DIRDELIMSTR, dbname);
  636.       strcpy(pathname, filename);
  637.    }
  639.    sprintf(filename, "%s.loo", pathname);
  640.    lut_mmap = Nlm_MemMapInit(filename);
  641.    if (lut_mmap == NULL)
  642.       ErrPostEx(SEV_FATAL, 1, 0, "Cannot map RPS BLAST lookup file");
  643.    info->lookup_header = (RPSLookupFileHeader *)lut_mmap->mmp_begin;
  645.    sprintf(filename, "%s.rps", pathname);
  646.    pssm_mmap = Nlm_MemMapInit(filename);
  647.    if (pssm_mmap == NULL)
  648.       ErrPostEx(SEV_FATAL, 1, 0, "Cannot map RPS BLAST profile file");
  649.    info->profile_header = (RPSProfileHeader *)pssm_mmap->mmp_begin;
  651.    num_db_seqs = info->profile_header->num_profiles;
  653.    sprintf(filename, "%s.aux", pathname);
  654.    auxfile = FileOpen(filename, "r");
  655.    if (auxfile == NULL)
  656.       ErrPostEx(SEV_FATAL, 1, 0,"Cannot open RPS BLAST parameters file");
  658.    fscanf(auxfile, "%s", buffer);
  659.    info->aux_info.orig_score_matrix = strdup(buffer);
  660.    fscanf(auxfile, "%d", &info->aux_info.gap_open_penalty);
  661.    fscanf(auxfile, "%d", &info->aux_info.gap_extend_penalty);
  662.    fscanf(auxfile, "%le", &info->aux_info.ungapped_k);
  663.    fscanf(auxfile, "%le", &info->aux_info.ungapped_h);
  664.    fscanf(auxfile, "%d", &info->aux_info.max_db_seq_length);
  665.    fscanf(auxfile, "%d", &info->aux_info.db_length);
  666.    fscanf(auxfile, "%lf", &info->aux_info.scale_factor);
  668.    info->aux_info.karlin_k = (double *)malloc(num_db_seqs * sizeof(double));
  669.    for (i = 0; i < num_db_seqs && !feof(auxfile); i++) {
  670.       fscanf(auxfile, "%d", &seq_size); /* not used */
  671.       fscanf(auxfile, "%le", &info->aux_info.karlin_k[i]);
  672.    }
  674.    if (i < num_db_seqs)
  675.       ErrPostEx(SEV_FATAL, 1, 0, "Missing Karlin parameters");
  677.    FileClose(auxfile);
  678.    *ppinfo = info;
  679.    *rps_mmap = lut_mmap;
  680.    *rps_pssm_mmap = pssm_mmap;
  681.    return 0;
  682. }
  684. int CBlastApplication::Run(void)
  685. {
  686.     Uint1 program_number;
  687.     int status = 0;
  689.     // Process command line args
  690.     const CArgs& args = GetArgs();
  691.     
  692.     BlastProgram2Number(args["program"].AsString().c_str(), &program_number);
  693.     EProgram program = static_cast<EProgram>(program_number);
  695.     Int4 strand_number = args["strand"].AsInteger();
  696.     ENa_strand strand;
  697.     Int4 from = args["qstart"].AsInteger();
  698.     Int4 to = args["qend"].AsInteger();
  699.     Int4 first_oid = 0;
  700.     Int4 last_oid = 0;
  702.     if (program == eBlastn ||
  703.         program == eBlastx ||
  704.         program == eRPSTblastn) {
  705.         if (strand_number == 1)
  706.             strand = eNa_strand_plus;
  707.         else if (strand_number == 2)
  708.             strand = eNa_strand_minus;
  709.         else
  710.             strand = eNa_strand_both;
  711.     }
  712.     else {
  713.         strand = eNa_strand_unknown;
  714.     }
  716.     InitScope();
  718.     int id_counter = 0;
  719.     // Read the query(ies) from input file; perform the setup
  720.     TSeqLocVector query_loc = 
  721.         BLASTGetSeqLocFromStream(args["query"].AsInputFile(),
  722.             *m_ObjMgr, strand, from, to, &id_counter,
  723.             args["lcase"].AsBoolean());
  725.     if (args["dbrange"]) {
  726.         const char* delimiters = " ,:;";
  727.         char* range_str = strdup(args["dbrange"].AsString().c_str());
  728.         first_oid = atoi(strtok(range_str, delimiters));
  729.         last_oid = atoi(strtok(NULL, delimiters));
  730.         sfree(range_str);
  731.     }
  732.     
  733.     bool db_is_na = (program == eBlastn || program == eTblastn || 
  734.                      program == eTblastx);
  736.     BlastSeqSrc* seq_src = 
  737.         SeqDbSrcInit(args["db"].AsString().c_str(), !db_is_na,
  738.                      first_oid, last_oid, NULL);
  740.     CBlastOptionsHandle* opts = CBlastOptionsFactory::Create(program);
  742.     Nlm_MemMapPtr rps_mmap = NULL;
  743.     Nlm_MemMapPtr rps_pssm_mmap = NULL;
  744.     RPSInfo *rps_info = NULL;
  745.     // Need to set up the RPS database information structure too
  746.     if (program == eRPSBlast || program == eRPSTblastn) {
  747.         if (BLAST_FillRPSInfo(&rps_info, &rps_mmap,
  748.             &rps_pssm_mmap, args["db"].AsString().c_str()) != 0) 
  749.         {
  750.             NCBI_THROW(CBlastException, eBadParameter, 
  751.                        "Cannot initialize RPS BLAST database");
  752.         }        
  753.     }
  754.     
  755.     ProcessCommandLineArgs(opts, seq_src, rps_info);
  757.     CDbBlast blaster(query_loc, seq_src, *opts, rps_info);
  759.     TSeqAlignVector seqalignv = blaster.Run();
  761.     BlastSeqSrcFree(seq_src);
  762.     if (rps_info) {
  763.         Nlm_MemMapFini(rps_mmap);
  764.         Nlm_MemMapFini(rps_pssm_mmap);
  765.         sfree(rps_info->aux_info.karlin_k);
  766.         sfree(rps_info->aux_info.orig_score_matrix);
  767.         sfree(rps_info);
  768.     }
  770.     FormatResults(blaster, seqalignv);
  772.     return status;
  773. }
  775. void CBlastApplication::Exit(void)
  776. {
  777.     SetDiagStream(0);
  778. }
  781. int main(int argc, const char* argv[] /*, const char* envp[]*/)
  782. {
  783.     return CBlastApplication().AppMain(argc, argv, 0, eDS_Default, 0);
  784. }