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

生物技术

开发平台:

C/C++

  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: omssacl.cpp,v $
  4.  * PRODUCTION Revision 1000.4  2004/06/01 18:09:13  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.13
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  9. /* $Id: omssacl.cpp,v 1000.4 2004/06/01 18:09:13 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 authors in any work or product based on this material.
  31.  *
  32.  * ===========================================================================
  33.  *
  34.  * Author:  Lewis Y. Geer
  35.  *  
  36.  * File Description:
  37.  *    command line OMSSA search
  38.  *
  39.  *
  40.  * ===========================================================================
  41.  */
  42. #include <ncbi_pch.hpp>
  43. #include <corelib/ncbiargs.hpp>
  44. #include <corelib/ncbiapp.hpp>
  45. #include <corelib/ncbienv.hpp>
  46. #include <corelib/ncbistre.hpp>
  47. #include <serial/serial.hpp>
  48. #include <serial/objistrasn.hpp>
  49. #include <serial/objistrasnb.hpp>
  50. #include <serial/objostrasn.hpp>
  51. #include <serial/objostrasnb.hpp>
  52. #include <serial/iterator.hpp>
  53. #include <objects/omssa/omssa__.hpp>
  54. #include <serial/objostrxml.hpp>
  55. #include <fstream>
  56. #include <string>
  57. #include <list>
  58. #include <stdio.h>
  59. #include "omssa.hpp"
  60. #include "SpectrumSet.hpp"
  61. #include "Mod.hpp"
  62. #include <readdb.h>
  63. USING_NCBI_SCOPE;
  64. USING_SCOPE(objects);
  65. USING_SCOPE(omssa);
  66. /////////////////////////////////////////////////////////////////////////////
  67. //
  68. //  COMSSA
  69. //
  70. //  Main application
  71. //
  72. typedef list <string> TStringList;
  73. class COMSSA : public CNcbiApplication {
  74. public:
  75.     virtual int Run();
  76.     virtual void Init();
  77.     void PrintMods(void);
  78.     void InsertMods(TStringList& List, CMSRequest::TVariable& Mods);
  79. };
  80. ///
  81. /// Insert modifications
  82. ///
  83. void COMSSA::InsertMods(TStringList& List, CMSRequest::TVariable& Mods)
  84. {
  85.     TStringList::iterator iList(List.begin());
  86.     int ModNum;
  87.     for(;iList != List.end(); iList++) {
  88. try {
  89.     ModNum = NStr::StringToInt(*iList);
  90. } catch (CStringException &e) {
  91.     ERR_POST(Info << "omssacl: unrecognized modification number");
  92.     continue;
  93. }
  94. Mods.push_back(ModNum);
  95.     }
  96. }
  97. ///
  98. /// print out a list of modification that can be used in OMSSA
  99. ///
  100. void COMSSA::PrintMods(void)
  101. {
  102.     int i;
  103.     for(i = 0; i < kNumMods; i++)
  104. cout << i << ": " << kModNames[i] << endl;
  105. }
  106. void COMSSA::Init()
  107. {
  108.     auto_ptr<CArgDescriptions> argDesc(new CArgDescriptions);
  109.     argDesc->AddDefaultKey("d", "blastdb", "Blast sequence library to search",
  110.    CArgDescriptions::eString, "nr");
  111.     argDesc->AddDefaultKey("f", "infile", "single dta file to search",
  112.    CArgDescriptions::eString, "");
  113.     argDesc->AddDefaultKey("fx", "xmlinfile", "multiple xml-encapsulated dta files to search",
  114.    CArgDescriptions::eString, "");
  115.     argDesc->AddDefaultKey("fb", "dtainfile", "multiple dta files separated by blank lines to search",
  116.    CArgDescriptions::eString, "");
  117.     argDesc->AddDefaultKey("o", "outfile", "asn.1 format search results filenam",
  118.    CArgDescriptions::eString, "");
  119.     argDesc->AddDefaultKey("ox", "outfile", "xml format search results filename",
  120.    CArgDescriptions::eString, "");
  121. #if MSGRAPH
  122.     argDesc->AddDefaultKey("g", "graph", "graph file to write out",
  123.    CArgDescriptions::eString, "");
  124. #endif
  125.     argDesc->AddDefaultKey("to", "pretol", "product ion mass tolerance in Da",
  126.    CArgDescriptions::eDouble, "0.8");
  127.     argDesc->AddDefaultKey("te", "protol", "precursor ion  mass tolerance in Da",
  128.    CArgDescriptions::eDouble, "2.0");
  129.     argDesc->AddDefaultKey("cl", "cutlo", "low intensity cutoff in % of max peak",
  130.    CArgDescriptions::eDouble, "0.0");
  131.     argDesc->AddDefaultKey("ch", "cuthi", "high intensity cutoff in % of max peak",
  132.    CArgDescriptions::eDouble, "0.2");
  133.     argDesc->AddDefaultKey("ci", "cutinc", "intensity cutoff increment in % of max peak",
  134.    CArgDescriptions::eDouble, "0.0005");
  135.     //    argDesc->AddDefaultKey("u", "cull", 
  136.     //    "number of peaks to leave in each 100Da bin",
  137.     //    CArgDescriptions::eInteger, "3");
  138.     argDesc->AddDefaultKey("v", "cleave", 
  139.    "number of missed cleavages allowed",
  140.    CArgDescriptions::eInteger, "1");
  141.     //    argDesc->AddDefaultKey("n", "numseq", 
  142.     //    "number of sequences to search (0 = all)",
  143.     //    CArgDescriptions::eInteger, "0");
  144.     argDesc->AddDefaultKey("x", "taxid", 
  145.    "taxid to search (0 = all)",
  146.    CArgDescriptions::eInteger, "0");
  147.     argDesc->AddDefaultKey("w1", "window1", 
  148.    "single charge window in Da",
  149.    CArgDescriptions::eInteger, "20");
  150.     argDesc->AddDefaultKey("w2", "window2", 
  151.    "double charge window in Da",
  152.    CArgDescriptions::eInteger, "14");
  153.     argDesc->AddDefaultKey("h1", "hit1", 
  154.    "number of peaks allowed in single charge window",
  155.    CArgDescriptions::eInteger, "2");
  156.     argDesc->AddDefaultKey("h2", "hit2", 
  157.    "number of peaks allowed in double charge window",
  158.    CArgDescriptions::eInteger, "2");
  159.     argDesc->AddDefaultKey("hl", "hitlist", 
  160.    "maximum number of hits retained for one spectrum",
  161.    CArgDescriptions::eInteger, "30");
  162.     argDesc->AddDefaultKey("mf", "fixedmod", 
  163.    "comma delimited list of id numbers for fixed modifications",
  164.    CArgDescriptions::eString,
  165.    "");
  166.     //  NStr::IntToString(eMSMod_Ccarbamidomethyl));
  167.     argDesc->AddDefaultKey("mv", "variablemod", 
  168.    "comma delimited list of id numbers for variable modifications",
  169.    CArgDescriptions::eString,
  170.    "");
  171.     // NStr::IntToString(eMSMod_Moxy));
  172.     argDesc->AddFlag("ml", "print a list of modifications and their corresponding id number");
  173.     SetupArgDescriptions(argDesc.release());
  174.     // allow info posts to be seen
  175.     SetDiagPostLevel(eDiag_Info);
  176. }
  177. int main(int argc, const char* argv[]) 
  178. {
  179.     COMSSA theTestApp;
  180.     return theTestApp.AppMain(argc, argv);
  181. }
  182. int COMSSA::Run()
  183. {    
  184.     try {
  185. CArgs args = GetArgs();
  186. // print out the modification list
  187. if(args["ml"]) {
  188.     PrintMods();
  189.     return 0;
  190. }
  191. _TRACE("omssa: initializing score");
  192. CSearch Search;
  193. int retval = Search.InitBlast(args["d"].AsString().c_str(), false);
  194. if(retval) {
  195.     ERR_POST(Fatal << "ommsacl: unable to initialize blastdb, error " 
  196.      << retval);
  197.     return 1;
  198. }
  199. _TRACE("ommsa: score initialized");
  200. CRef <CSpectrumSet> Spectrumset(new CSpectrumSet);
  201. if(args["fx"].AsString().size() != 0) {
  202.     ifstream PeakFile(args["fx"].AsString().c_str());
  203.     if(!PeakFile) {
  204. ERR_POST(Fatal <<" omssacl: not able to open spectrum file " <<
  205.  args["fx"].AsString());
  206. return 1;
  207.     }
  208.     if(Spectrumset->LoadMultDTA(PeakFile) != 0) {
  209. ERR_POST(Fatal <<" omssacl: error reading spectrum file " <<
  210.  args["fx"].AsString());
  211. return 1;
  212.     }
  213. }
  214. else if(args["f"].AsString().size() != 0) {
  215.     ifstream PeakFile(args["f"].AsString().c_str());
  216.     if(!PeakFile) {
  217. ERR_POST(Fatal << "omssacl: not able to open spectrum file " <<
  218.  args["f"].AsString());
  219. return 1;
  220.     }
  221.     if(Spectrumset->LoadDTA(PeakFile) != 0) {
  222. ERR_POST(Fatal << "omssacl: not able to read spectrum file " <<
  223.  args["f"].AsString());
  224. return 1;
  225.     }
  226. }
  227. else if(args["fb"].AsString().size() != 0) {
  228.     ifstream PeakFile(args["fb"].AsString().c_str());
  229.     if(!PeakFile) {
  230. ERR_POST(Fatal << "omssacl: not able to open spectrum file " <<
  231.  args["fb"].AsString());
  232. return 1;
  233.     }
  234.     if(Spectrumset->LoadMultBlankLineDTA(PeakFile) != 0) {
  235. ERR_POST(Fatal << "omssacl: not able to read spectrum file " <<
  236.  args["fb"].AsString());
  237. return 1;
  238.     }
  239. }
  240. else {
  241.     ERR_POST(Fatal << "omssatest: input file not given.");
  242.     return 1;
  243. }
  244. CMSResponse Response;
  245. CMSRequest Request;
  246. Request.SetSearchtype(eMSSearchType_monoisotopic);
  247. Request.SetPeptol(args["te"].AsDouble());
  248. Request.SetMsmstol(args["to"].AsDouble());
  249. Request.SetCutlo(args["cl"].AsDouble());
  250. Request.SetCuthi(args["ch"].AsDouble());
  251. Request.SetCutinc(args["ci"].AsDouble());
  252. Request.SetSinglewin(args["w1"].AsInteger());
  253. Request.SetDoublewin(args["w2"].AsInteger());
  254. Request.SetSinglenum(args["h1"].AsInteger());
  255. Request.SetDoublenum(args["h2"].AsInteger());
  256. Request.SetEnzyme(eMSEnzymes_trypsin);
  257. Request.SetMissedcleave(args["v"].AsInteger());
  258. Request.SetSpectra(*Spectrumset);
  259. {
  260.     TStringList List;
  261.     NStr::Split(args["mv"].AsString(), ",", List);
  262.     InsertMods(List, Request.SetVariable());
  263.     List.clear();
  264.     NStr::Split(args["mf"].AsString(), ",", List);
  265.     InsertMods(List, Request.SetFixed());
  266. }
  267. Request.SetDb(args["d"].AsString());
  268. // Request.SetCull(args["u"].AsInteger());
  269. Request.SetHitlistlen(args["hl"].AsInteger());
  270. if(args["x"].AsInteger() != 0) 
  271.     Request.SetTaxids().push_back(args["x"].AsInteger());
  272. _TRACE("omssa: search begin");
  273. Search.Search(Request, Response);
  274. _TRACE("omssa: search end");
  275. #if _DEBUG
  276. // read out hits
  277. CMSResponse::THitsets::const_iterator iHits;
  278. iHits = Response.GetHitsets().begin();
  279. for(; iHits != Response.GetHitsets().end(); iHits++) {
  280.     CRef< CMSHitSet > HitSet =  *iHits;
  281.     ERR_POST(Info << "Hitset: " << HitSet->GetNumber());
  282.     if( HitSet-> CanGetError() && HitSet->GetError() ==
  283. eMSHitError_notenuffpeaks) {
  284. ERR_POST(Info << "Hitset Empty");
  285. continue;
  286.     }
  287.     CRef< CMSHits > Hit;
  288.     CMSHitSet::THits::const_iterator iHit; 
  289.     CMSHits::TPephits::const_iterator iPephit;
  290.     for(iHit = HitSet->GetHits().begin();
  291. iHit != HitSet->GetHits().end(); iHit++) {
  292. ERR_POST(Info << (*iHit)->GetPepstring() << ": " << "P = " << 
  293.  (*iHit)->GetPvalue() << " E = " <<
  294.  (*iHit)->GetEvalue());    
  295. for(iPephit = (*iHit)->GetPephits().begin();
  296.     iPephit != (*iHit)->GetPephits().end();
  297.     iPephit++) {
  298.     ERR_POST(Info << (*iPephit)->GetGi() << ": " << (*iPephit)->GetStart() << "-" << (*iPephit)->GetStop() << ":" << (*iPephit)->GetDefline());
  299. }
  300.     
  301.     }
  302. }
  303. #endif
  304. if(args["o"].AsString() != "") {
  305.     auto_ptr<CObjectOStream>
  306. txt_out(CObjectOStream::Open(args["o"].AsString().c_str(),
  307.      eSerial_AsnText));
  308.     txt_out->Write(ObjectInfo(Response));
  309. }
  310. else if(args["ox"].AsString() != "") {
  311.     CNcbiOfstream os(args["ox"].AsString().c_str());
  312.     auto_ptr<CObjectOStreamXml>
  313. txt_out(new CObjectOStreamXml(os, false));
  314.     //txt_out->SetEnforcedStdXml();
  315.     txt_out->Write(ObjectInfo(Response));
  316. }
  317. #if MSGRAPH
  318. if(args["g"].AsString() != "") {
  319.     CMSDrawSpectra DrawIt(kMSWIDTH, kMSHEIGHT);
  320.     DrawIt.Init(Request);
  321.     DrawIt.DrawHit(Response, 0, Search, true, true);
  322.     DrawIt.DrawSpectra();
  323.     FILE *os;
  324.     os = fopen(args["g"].AsString().c_str(), "wb");
  325.     DrawIt.Output(os);
  326.     fclose(os);
  327. }
  328. #endif
  329.     } catch (NCBI_NS_STD::exception& e) {
  330. ERR_POST(Fatal << "Exception in COMSSA::Run: " << e.what());
  331.     }
  332.     return 0;
  333. }
  334. /*
  335.   $Log: omssacl.cpp,v $
  336.   Revision 1000.4  2004/06/01 18:09:13  gouriano
  337.   PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.13
  338.   Revision 1.13  2004/05/27 20:52:15  lewisg
  339.   better exception checking, use of AutoPtr, command line parsing
  340.   Revision 1.12  2004/05/21 21:41:03  gorelenk
  341.   Added PCH ncbi_pch.hpp
  342.   Revision 1.11  2004/03/30 19:36:59  lewisg
  343.   multiple mod code
  344.   Revision 1.10  2004/03/16 22:09:11  gorelenk
  345.   Changed include for private header.
  346.   Revision 1.9  2004/03/01 18:24:08  lewisg
  347.   better mod handling
  348.   Revision 1.8  2003/12/04 23:39:09  lewisg
  349.   no-overlap hits and various bugfixes
  350.   Revision 1.7  2003/11/20 15:40:53  lewisg
  351.   fix hitlist bug, change defaults
  352.   Revision 1.6  2003/11/18 18:16:04  lewisg
  353.   perf enhancements, ROCn adjusted params made default
  354.   Revision 1.5  2003/11/13 19:07:38  lewisg
  355.   bugs: iterate completely over nr, don't initialize blastdb by iteration
  356.   Revision 1.4  2003/11/10 22:24:12  lewisg
  357.   allow hitlist size to vary
  358.   Revision 1.3  2003/10/27 20:10:55  lewisg
  359.   demo program to read out omssa results
  360.   Revision 1.2  2003/10/21 21:12:17  lewisg
  361.   reorder headers
  362.   Revision 1.1  2003/10/20 21:32:13  lewisg
  363.   ommsa toolkit version
  364.   Revision 1.1  2003/10/07 18:02:28  lewisg
  365.   prep for toolkit
  366.   Revision 1.10  2003/10/06 18:14:17  lewisg
  367.   threshold vary
  368.   Revision 1.9  2003/07/21 20:25:03  lewisg
  369.   fix missing peak bug
  370.   Revision 1.8  2003/07/17 18:45:50  lewisg
  371.   multi dta support
  372.   Revision 1.7  2003/07/07 16:17:51  lewisg
  373.   new poisson distribution and turn off histogram
  374.   Revision 1.6  2003/05/01 14:52:10  lewisg
  375.   fixes to scoring
  376.   Revision 1.5  2003/04/18 20:46:52  lewisg
  377.   add graphing to omssa
  378.   Revision 1.4  2003/04/04 18:43:51  lewisg
  379.   tax cut
  380.   Revision 1.3  2003/04/02 18:49:51  lewisg
  381.   improved score, architectural changes
  382.   Revision 1.2  2003/02/10 19:37:56  lewisg
  383.   perf and web page cleanup
  384.   Revision 1.1  2003/02/03 20:39:06  lewisg
  385.   omssa cgi
  386. */