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

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

生物技术

开发平台：

C/C++

omssa.cpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: omssa.cpp,v $
  4.  * PRODUCTION Revision 1000.4  2004/06/01 18:09:10  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.20
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /* $Id: omssa.cpp,v 1000.4 2004/06/01 18:09:10 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 authors in any work or product based on this material.
  32. *
  33. * ===========================================================================
  34. *
  35. * Author:  Lewis Y. Geer
  36. *  
  37. * File Description:
  38. *    code to do the ms/ms search and score matches
  39. *
  40. * ===========================================================================
  41. */
  43. #include <ncbi_pch.hpp>
  44. #include <corelib/ncbistd.hpp>
  45. #include <corelib/ncbiargs.hpp>
  46. #include <corelib/ncbiapp.hpp>
  47. #include <corelib/ncbienv.hpp>
  48. #include <corelib/ncbidiag.hpp>
  49. #include <corelib/ncbi_limits.hpp>
  51. #include <fstream>
  52. #include <string>
  53. #include <list>
  54. #include <deque>
  55. #include <algorithm>
  56. #include <math.h>
  58. #include "SpectrumSet.hpp"
  59. #include "omssa.hpp"
  61. #include "nrutil.h"
  62. #include <ncbimath.h>
  64. USING_NCBI_SCOPE;
  65. USING_SCOPE(objects);
  66. USING_SCOPE(omssa);
  68. /////////////////////////////////////////////////////////////////////////////
  69. //
  70. //  CSearch::
  71. //
  72. //  Performs the ms/ms search
  73. //
  76. CSearch::CSearch(void): rdfp(0)
  77. {
  78. }
  80. int CSearch::InitBlast(const char *blastdb, bool InitDb)
  81. {
  82.     if(!blastdb) return 0;
  83.     if(rdfp) readdb_destruct(rdfp);
  84.     rdfp = readdb_new_ex2((char *)blastdb,TRUE, READDB_NEW_DO_TAXDB | 
  85.   READDB_NEW_INDEX, NULL, NULL);
  86.     if(!rdfp) return 1;
  88.     numseq = readdb_get_num_entries(rdfp);
  90.     if(InitDb) {
  91. int iSearch, length;
  92. char testchar;
  93. unsigned char *Sequence;
  95. for(iSearch = 0; iSearch < numseq; iSearch++) {
  96.     length = readdb_get_sequence(rdfp, iSearch, &Sequence);
  97.     testchar = Sequence[0];
  98. }
  99.     }
  100.     return 0;
  102. }
  105. // create the ladders from sequence
  107. int CSearch::CreateLadders(unsigned char *Sequence, int iSearch, int position,
  108.    int endposition,
  109.    int *Masses, int iMissed, CAA& AA, 
  110.    CLadder& BLadder,
  111.    CLadder& YLadder, CLadder& B2Ladder,
  112.    CLadder& Y2Ladder,
  113.    unsigned ModMask,
  114.    const char **Site,
  115.    int *DeltaMass,
  116.    int NumMod)
  117. {
  118.     if(!BLadder.CreateLadder(kBIon, 1, (char *)Sequence, iSearch,
  119.      position, endposition, Masses[iMissed], 
  120.      MassArray, AA, ModMask, Site, DeltaMass,
  121.      NumMod)) return 1;
  122.     if(!YLadder.CreateLadder(kYIon, 1, (char *)Sequence, iSearch,
  123.  position, endposition, Masses[iMissed], 
  124.      MassArray, AA, ModMask, Site, DeltaMass,
  125.      NumMod)) return 1;
  126.     B2Ladder.CreateLadder(kBIon, 2, (char *)Sequence, iSearch,
  127.   position, endposition, 
  128.   Masses[iMissed], 
  129.   MassArray, AA, ModMask, Site, DeltaMass,
  130.   NumMod);
  131.     Y2Ladder.CreateLadder(kYIon, 2, (char *)Sequence, iSearch,
  132.   position, endposition,
  133.   Masses[iMissed], 
  134.   MassArray, AA, ModMask, Site, DeltaMass,
  135.   NumMod);
  136.     
  137.     return 0;
  138. }
  141. // compare ladders to experiment
  142. int CSearch::CompareLadders(CLadder& BLadder,
  143.     CLadder& YLadder, CLadder& B2Ladder,
  144.     CLadder& Y2Ladder, CMSPeak *Peaks,
  145.     bool OrLadders,  TMassPeak *MassPeak)
  146. {
  147.     if(MassPeak && MassPeak->Charge >= kConsiderMult) {
  148. Peaks->CompareSorted(BLadder, MSCULLED2, 0); 
  149. Peaks->CompareSorted(YLadder, MSCULLED2, 0); 
  150. Peaks->CompareSorted(B2Ladder, MSCULLED2, 0); 
  151. Peaks->CompareSorted(Y2Ladder, MSCULLED2, 0);
  152. if(OrLadders) {
  153.     BLadder.Or(B2Ladder);
  154.     YLadder.Or(Y2Ladder);
  155. }
  156.     }
  157.     else {
  158. Peaks->CompareSorted(BLadder, MSCULLED1, 0); 
  159. Peaks->CompareSorted(YLadder, MSCULLED1, 0); 
  160.     }
  161.     return 0;
  162. }
  164. // compare ladders to experiment
  165. bool CSearch::CompareLaddersTop(CLadder& BLadder,
  166.     CLadder& YLadder, CLadder& B2Ladder,
  167.     CLadder& Y2Ladder, CMSPeak *Peaks,
  168.     TMassPeak *MassPeak)
  169. {
  170.     if(MassPeak && MassPeak->Charge >=  kConsiderMult ) {
  171. if(Peaks->CompareTop(BLadder)) return true; 
  172. if(Peaks->CompareTop(YLadder)) return true; 
  173. if(Peaks->CompareTop(B2Ladder)) return true; 
  174. if(Peaks->CompareTop(Y2Ladder)) return true;
  175.     }
  176.     else {
  177. if(Peaks->CompareTop(BLadder)) return true; 
  178. if(Peaks->CompareTop(YLadder)) return true; 
  179.     }
  180.     return false;
  181. }
  183. #ifdef _DEBUG
  184. #define CHECKGI
  185. #ifdef CHECKGI
  186. bool CheckGi(int gi)
  187. {
  188.     if(gi == 21742354) {
  189. ERR_POST(Info << "test seq");
  190. return true;
  191.     }
  192.     return false;
  193. }
  194. #endif
  195. #endif
  197. // loads spectra into peaks
  198. void CSearch::Spectrum2Peak(CMSRequest& MyRequest, CMSPeakSet& PeakSet)
  199. {
  200.     CSpectrumSet::Tdata::const_iterator iSpectrum;
  201.     CMSPeak* Peaks;
  203.     iSpectrum = MyRequest.GetSpectra().Get().begin();
  204.     for(; iSpectrum != MyRequest.GetSpectra().Get().end(); iSpectrum++) {
  205. CRef <CMSSpectrum> Spectrum =  *iSpectrum;
  206. if(!Spectrum) {
  207.     ERR_POST(Error << "omssa: unable to find spectrum");
  208.     return;
  209. }
  211. Peaks = new CMSPeak(MyRequest.GetHitlistlen());
  212. if(!Peaks) {
  213.     ERR_POST(Error << "omssa: unable to allocate CMSPeak");
  214.     return;
  215. }
  216. if(Peaks->Read(*Spectrum, MyRequest.GetMsmstol()) != 0) {
  217.     ERR_POST(Error << "omssa: unable to read spectrum into CMSPeak");
  218.     return;
  219. }
  220. if(Spectrum->CanGetName()) Peaks->SetName() = Spectrum->GetName();
  221. if(Spectrum->CanGetNumber())
  222.     Peaks->SetNumber() = Spectrum->GetNumber();
  224. Peaks->Sort();
  225. Peaks->SetComputedCharge(3);
  226. Peaks->InitHitList();
  227. Peaks->CullAll(MyRequest.GetCutlo(), MyRequest.GetSinglewin(),
  228.        MyRequest.GetDoublewin(), MyRequest.GetSinglenum(),
  229.        MyRequest.GetDoublenum());
  231. if(Peaks->GetNum(MSCULLED1) < MSPEAKMIN)  {
  232.     ERR_POST(Info << "omssa: not enough peaks in spectra");
  233.     Peaks->SetError(eMSHitError_notenuffpeaks);
  234. }
  235. PeakSet.AddPeak(Peaks);
  237.     }
  238.     PeakSet.SortPeaks(static_cast <int> (MyRequest.GetPeptol()*MSSCALE));
  240. }
  242. // compares TMassMasks.  Lower m/z first in sort.
  243. struct CMassMaskCompare {
  244.     bool operator() (const TMassMask& x, const TMassMask& y)
  245.     {
  246. if(x.Mass < y.Mass) return true;
  247. return false;
  248.     }
  249. };
  251. // update sites and masses for new peptide
  252. void CSearch::UpdateWithNewPep(int Missed,
  253.        const char *PepStart[],
  254.        const char *PepEnd[], 
  255.        int NumMod[], 
  256.        const char *Site[][MAXMOD],
  257.        int DeltaMass[][MAXMOD],
  258.        int Masses[],
  259.        int EndMasses[])
  260. {
  261.     // iterate over missed cleavages
  262.     int iMissed;
  263.     // maximum mods allowed
  264.     int ModMax; 
  265.     // iterate over mods
  266.     int iMod;
  267.     
  269.     // update the longer peptides to add the new peptide (Missed-1) on the end
  270.     for(iMissed = 0; iMissed < Missed - 1; iMissed++) {
  271. // skip start
  272. if(PepStart[iMissed] == (const char *)-1) continue;
  273. // reset the end sequences
  274. PepEnd[iMissed] = PepEnd[Missed - 1];
  276. // update new mod masses to add in any new mods from new peptide
  278. // first determine the maximum value for updated mod list
  279. if(NumMod[iMissed] + NumMod[Missed-1] >= MAXMOD)
  280.     ModMax = MAXMOD - NumMod[iMissed];
  281. else ModMax = NumMod[Missed-1];
  283. // now interatu thru the new entries
  284. for(iMod = NumMod[iMissed]; 
  285.     iMod < NumMod[iMissed] + ModMax;
  286.     iMod++) {
  287.     Site[iMissed][iMod] = 
  288. Site[Missed-1][iMod - NumMod[iMissed]];
  290.     DeltaMass[iMissed][iMod] = 
  291. DeltaMass[Missed-1][iMod - NumMod[iMissed]];
  293.     //     Masses[iMissed][iMod] = 
  294.     // Masses[Missed-1][iMod - NumMod[iMissed] + 1] +
  295.     // Masses[iMissed][NumMod[iMissed] - 1];
  296. }
  298. // update old masses
  299. Masses[iMissed] += Masses[Missed - 1];
  301. // update end masses
  302. EndMasses[iMissed] = EndMasses[Missed - 1];
  304. // update number of Mods
  305. NumMod[iMissed] += ModMax;
  306.     }     
  307. }
  310. // create the various combinations of mods
  311. void CSearch::CreateModCombinations(int Missed,
  312.     const char *PepStart[],
  313.     TMassMask MassAndMask[][MAXMOD2],
  314.     int Masses[],
  315.     int EndMasses[],
  316.     int NumMod[],
  317.     int DeltaMass[][MAXMOD],
  318.     unsigned NumMassAndMask[])
  319. {
  320.     // need to iterate thru combinations that have iMod.
  321.     // i.e. iMod = 3 and NumMod=5
  322.     // 00111, 01011, 10011, 10101, 11001, 11010, 11100, 01101,
  323.     // 01110
  324.     // i[0] = 0 --> 5-3, i[1] = i[0]+1 -> 5-2, i[3] = i[1]+1 -> 5-1
  325.     // then construct bool mask
  326.     
  327.     // holders for calculated modification mask and modified peptide masses
  328.     unsigned Mask, MassOfMask;
  329.     // iterate thru active mods
  330.     int iiMod;
  331.     // keep track of the number of unique masks created.  each corresponds to a ladder
  332.     int iModCount;
  333.     // missed cleavage
  334.     int iMissed;
  335.     // number of mods to consider
  336.     int iMod;
  337.     // positions of mods
  338.     int ModIndex[MAXMOD];
  340.     // go thru missed cleaves
  341.     for(iMissed = 0; iMissed < Missed; iMissed++) {
  342. // skip start
  343. if(PepStart[iMissed] == (const char *)-1) continue;
  344. iModCount = 0;
  346. // set up non-modified mass
  347. MassAndMask[iMissed][iModCount].Mass = 
  348.     Masses[iMissed] + EndMasses[iMissed];
  349. MassAndMask[iMissed][iModCount].Mask = 0;
  350. iModCount++;
  352. // go thru number of mods allowed
  353. for(iMod = 0; iMod < NumMod[iMissed] && iModCount < MAXMOD2; iMod++) {
  354.     
  355.     // initialize ModIndex that points to mod sites
  356.     InitModIndex(ModIndex, iMod);
  357.     do {
  359. // calculate mass
  360. MassOfMask = Masses[iMissed] + EndMasses[iMissed];
  361. for(iiMod = 0; iiMod <= iMod; iiMod++ ) 
  362.     MassOfMask += DeltaMass[iMissed][ModIndex[iiMod]];
  363. // make bool mask
  364. Mask = MakeBoolMask(ModIndex, iMod);
  365. // put mass and mask into storage
  366. MassAndMask[iMissed][iModCount].Mass = MassOfMask;
  367. MassAndMask[iMissed][iModCount].Mask = Mask;
  368. // keep track of the  number of ladders
  369. iModCount++;
  371.     } while(iModCount < MAXMOD2 &&
  372.     CalcModIndex(ModIndex, iMod, NumMod[iMissed]));
  373. } // iMod
  375. // sort mask and mass by mass
  376. sort(MassAndMask[iMissed], MassAndMask[iMissed] + iModCount,
  377.      CMassMaskCompare());
  378. // keep track of number of MassAndMask
  379. NumMassAndMask[iMissed] = iModCount;
  381.     } // iMissed
  382. }
  384. //#define MSSTATRUN
  386. int CSearch::Search(CMSRequest& MyRequest, CMSResponse& MyResponse)
  387. {
  389.     try {
  391. int length;
  392. CTrypsin trypsin;
  393. unsigned header;
  395. #ifdef CHECKGI
  396. char *blastdefline;
  397. SeqId *sip;
  398. #endif
  400. CLadder BLadder[MAXMOD2], YLadder[MAXMOD2], B2Ladder[MAXMOD2],
  401.     Y2Ladder[MAXMOD2];
  402. bool LadderCalc[MAXMOD2];  // have the ladders been calculated?
  403. CAA AA;
  405. int Missed = MyRequest.GetMissedcleave()+1;  // number of missed cleaves allowed + 1
  406. int iMissed; // iterate thru missed cleavages
  407.     
  408. int iSearch, hits;
  409. unsigned char *Sequence;  // start position
  410. int endposition, position;
  411. MassArray.Init(MyRequest.GetFixed(), MyRequest.GetSearchtype());
  412. VariableMods.Init(MyRequest.GetVariable());
  413. const int *IntMassArray = MassArray.GetIntMass();
  414. const char *PepStart[MAXMISSEDCLEAVE];
  415. const char *PepEnd[MAXMISSEDCLEAVE];
  417. // the position within the peptide of a variable modification
  418. const char *Site[MAXMISSEDCLEAVE][MAXMOD];
  419. // the modification mass at the Site
  420. int DeltaMass[MAXMISSEDCLEAVE][MAXMOD];
  421. // the number of modified sites + 1 unmodified
  422. int NumMod[MAXMISSEDCLEAVE];
  425. // calculated masses and masks
  426. TMassMask MassAndMask[MAXMISSEDCLEAVE][MAXMOD2];
  427. // the number of masses and masks for each peptide
  428. unsigned NumMassAndMask[MAXMISSEDCLEAVE];
  430. // set up mass array, indexed by missed cleavage
  431. // note that EndMasses is the end mass of peptide, kept separate to allow
  432. // reuse of Masses array in missed cleavage calc
  433. int Masses[MAXMISSEDCLEAVE];
  434. int EndMasses[MAXMISSEDCLEAVE];
  437. int iMod;   // used to iterate thru modifications
  439. bool SequenceDone;  // are we done iterating through the sequences?
  440. int taxid;
  441. const CMSRequest::TTaxids& Tax = MyRequest.GetTaxids();
  442. CMSRequest::TTaxids::const_iterator iTax;
  443. CMSHit NewHit;  // a new hit of a ladder to an m/z value
  444. CMSHit *NewHitOut;  // copy of new hit
  447. CMSPeakSet PeakSet;
  448. TMassPeak *MassPeak; // peak currently in consideration
  449. CMSPeak* Peaks;
  451. Spectrum2Peak(MyRequest, PeakSet);
  453. #ifdef MSSTATRUN
  454. ofstream charge1("charge1.txt");
  455. ofstream charge2("charge2.txt");
  456. ofstream charge3("charge3.txt");
  457. #endif
  459. // iterate through sequences
  460. for(iSearch = 0; iSearch < numseq; iSearch++) {
  461.     if(iSearch/10000*10000 == iSearch) ERR_POST(Info << "sequence " << 
  462. iSearch);
  463.     header = 0;
  464. #ifdef CHECKGI
  465.     int testgi;
  466. #endif 
  467.     if(MyRequest.IsSetTaxids()) {
  468. while(readdb_get_header_ex(rdfp, iSearch, &header, 
  469. #ifdef CHECKGI
  470.    &sip, &blastdefline,
  471. #else
  472.    NULL, NULL,
  473. #endif
  475.    &taxid, NULL, NULL)) {
  477. #ifdef CHECKGI
  478.     SeqId *bestid = SeqIdFindBest(sip, SEQID_GI);
  479.     CheckGi(bestid->data.intvalue);
  480.     testgi = bestid->data.intvalue;
  481.     //     cout << testgi << endl;
  482. #endif 
  484. #ifdef CHECKGI
  485.     MemFree(blastdefline);
  486.     SeqIdSetFree(sip);
  487. #endif
  488.     for(iTax = Tax.begin(); iTax != Tax.end(); iTax++) {
  489. if(taxid == *iTax) goto TaxContinue;
  490.     } 
  491.     //else goto TaxContinue;
  492. }
  493. continue;
  494.     }
  495. TaxContinue:
  496.     length = readdb_get_sequence(rdfp, iSearch, &Sequence);
  497.     SequenceDone = false;
  498.     // PepStart = (const char *)Sequence[0];
  500.     // initialize missed cleavage matrix
  501.     for(iMissed = 0; iMissed < Missed; iMissed++) {
  502. PepStart[iMissed] = (const char *)-1; // mark start
  503. PepEnd[iMissed] = (const char *)Sequence;
  504. Masses[iMissed] = 0;
  505. EndMasses[iMissed] = 0;
  506. NumMod[iMissed] = 0;
  508. DeltaMass[iMissed][0] = 0;
  509. Site[iMissed][0] = (const char *)-1;
  510.     }
  511.     PepStart[Missed - 1] = (const char *)Sequence;
  513.     // iterate thru the sequence by digesting it
  514.     while(!SequenceDone) {
  516. // zero out no missed cleavage peptide mass and mods
  517. // note that Masses and EndMass are separate to reuse
  518. // masses during the missed cleavage calculation
  519. Masses[Missed - 1] = 0;
  520. EndMasses[Missed - 1] = 0;
  521. NumMod[Missed - 1] = 0;
  522. // init no modification elements
  523. Site[Missed - 1][0] = (const char *)-1;
  524. DeltaMass[Missed - 1][0] = 0;
  526. // calculate new stop and mass
  527. SequenceDone = 
  528.     trypsin.CalcAndCut((const char *)Sequence,
  529.        (const char *)Sequence + length - 1, 
  530.        &(PepEnd[Missed - 1]),
  531.        &(Masses[Missed - 1]),
  532.        NumMod[Missed - 1],
  533.        MAXMOD,
  534.        &(EndMasses[Missed - 1]),
  535.        VariableMods,
  536.        Site[Missed - 1],
  537.        DeltaMass[Missed - 1],
  538.        IntMassArray);
  541. UpdateWithNewPep(Missed, PepStart, PepEnd, NumMod, Site,
  542.  DeltaMass, Masses, EndMasses);
  544. CreateModCombinations(Missed, PepStart, MassAndMask, Masses,
  545.       EndMasses, NumMod, DeltaMass, NumMassAndMask);
  548. int OldMass;  // keeps the old peptide mass for comparison
  549. bool NoMassMatch;  // was there a match to the old mass?
  551. for(iMissed = 0; iMissed < Missed; iMissed++) {     
  552.     if(PepStart[iMissed] == (const char *)-1) continue;  // skip start
  554.     // get the start and stop position, inclusive, of the peptide
  555.     position =  PepStart[iMissed] - (const char *)Sequence;
  556.     endposition = PepEnd[iMissed] - (const char *)Sequence;
  558.     // init bool for "Has ladder been calculated?"
  559.     for(iMod = 0; iMod < MAXMOD2; iMod++) LadderCalc[iMod] = false;
  561.     OldMass = 0;
  562.     NoMassMatch = true;
  564.     // go thru total number of mods
  565.     for(iMod = 0; iMod < NumMassAndMask[iMissed]; iMod++) {
  566.     
  567. // have we seen this mass before?
  568. if(MassAndMask[iMissed][iMod].Mass == OldMass &&
  569.    NoMassMatch) continue;
  570. NoMassMatch = true;
  571. OldMass = MassAndMask[iMissed][iMod].Mass;
  572. // return peak where theoretical mass is < precursor mass + tol
  573. MassPeak = PeakSet.GetIndexLo(OldMass);
  575. for(;MassPeak < PeakSet.GetEndMassPeak() && 
  576. OldMass > MassPeak->Mass - MassPeak->Peptol;
  577.     MassPeak++) {
  578.     Peaks = MassPeak->Peak;
  579.     // make sure we look thru other mod masks with the same mass
  580.     NoMassMatch = false;
  582. #ifdef CHECKGI
  583.     SeqId *bestid;
  585.     header = 0;
  586.     readdb_get_header(rdfp, iSearch, &header, &sip, &blastdefline);
  587.     bestid = SeqIdFindBest(sip, SEQID_GI);
  588.     if(/*Peaks->GetNumber() == 245 && */CheckGi(bestid->data.intvalue))
  589. cerr << "hello" << endl;
  590.     // int testgi = bestid->data.intvalue;
  591.     MemFree(blastdefline);
  592.     SeqIdSetFree(sip);
  593. #endif
  596.     if(!LadderCalc[iMod]) {
  597. LadderCalc[iMod] = true;
  598. if(CreateLadders(Sequence, iSearch, position,
  599.  endposition,
  600.  Masses,
  601.  iMissed, AA,
  602.  BLadder[iMod],
  603.  YLadder[iMod],
  604.  B2Ladder[iMod],
  605.  Y2Ladder[iMod],
  606.  MassAndMask[iMissed][iMod].Mask,
  607.  Site[iMissed],
  608.  DeltaMass[iMissed],
  609.  NumMod[iMissed]) != 
  610.    0) continue;
  611. // continue to next sequence if ladders not successfully made
  612.     }
  613.     else {
  614. BLadder[iMod].ClearHits();
  615. YLadder[iMod].ClearHits();
  616. B2Ladder[iMod].ClearHits();
  617. Y2Ladder[iMod].ClearHits();
  618.     }
  620.     if(
  621. #ifdef MSSTATRUN
  622.        true
  623. #else
  624.        CompareLaddersTop(BLadder[iMod], 
  625.  YLadder[iMod],
  626.  B2Ladder[iMod], 
  627.  Y2Ladder[iMod],
  628.  Peaks, MassPeak)
  629. #endif
  630.        ) {
  631. // end of new addition
  634. Peaks->SetPeptidesExamined(MassPeak->Charge)++;
  635. #ifdef CHECKGI
  636. /*if(Peaks->GetNumber() == 245)*/
  637. //     cout << testgi << " " << position << " " << endposition << endl;
  638. #endif
  639. Peaks->ClearUsedAll();
  640. CompareLadders(BLadder[iMod], 
  641.        YLadder[iMod],
  642.        B2Ladder[iMod], 
  643.        Y2Ladder[iMod],
  644.        Peaks, false, MassPeak);
  645. hits = BLadder[iMod].HitCount() + 
  646.     YLadder[iMod].HitCount() +
  647.     B2Ladder[iMod].HitCount() +
  648.     Y2Ladder[iMod].HitCount();
  649. #ifdef MSSTATRUN
  650. switch (MassPeak->Charge)
  651.     {
  652.     case 1:
  653. charge1 << hits << endl;
  654. break;
  655.     case 2:
  656. charge2 << hits << endl;
  657. break;
  658.     case 3:
  659. charge3 << hits << endl;
  660. break;
  661.     default:
  662. break;
  663.     }
  664. #endif
  665. if(hits >= MSHITMIN) {
  666.     // need to save mods.  bool map?
  667.     NewHit.SetHits(hits);
  668.     NewHit.SetCharge(MassPeak->Charge);
  669.     // only record if hit kept
  670.     if(Peaks->AddHit(NewHit, NewHitOut)) {
  671. NewHitOut->SetStart(position);
  672. NewHitOut->SetStop(endposition);
  673. NewHitOut->SetSeqIndex(iSearch);
  674. NewHitOut->SetMass(MassPeak->Mass);
  675. // record the hits
  676. Peaks->ClearUsedAll();
  677. NewHitOut->
  678.     RecordMatches(BLadder[iMod],
  679.   YLadder[iMod],
  680.   B2Ladder[iMod], 
  681.   Y2Ladder[iMod],
  682.   Peaks);
  683.     }
  684. }
  685.     } // new addition
  686. } // MassPeak
  687.     } //iMod
  688. } // iMissed
  689. if(!SequenceDone) {
  690.     // get rid of longest peptide and move the other peptides down the line
  691.     for(iMissed = 0; iMissed < Missed - 1; iMissed++) {
  692. // move masses to next missed cleavage
  693. NumMod[iMissed] = NumMod[iMissed + 1];
  694. Masses[iMissed] = Masses[iMissed + 1];
  695. // don't move EndMasses as they are recalculated
  697. // move the modification data
  698. for(iMod = 0; iMod < NumMod[iMissed]; iMod++) {
  699.     DeltaMass[iMissed][iMod] = 
  700. DeltaMass[iMissed + 1][iMod];
  701.     Site[iMissed][iMod] = 
  702. Site[iMissed + 1][iMod];
  703. }
  705. // copy starts to next missed cleavage
  706. PepStart[iMissed] = PepStart[iMissed + 1];
  707.     }
  710.     // init new start from old stop
  711.     PepEnd[Missed-1] += 1;
  712.     PepStart[Missed-1] = PepEnd[Missed-1];
  713.     // PepStart = PepEnd + 1;
  714. } 
  715.     }
  716. }
  718. // read out hits
  719. SetResult(PeakSet, MyResponse, MyRequest.GetCutlo(), 
  720.   MyRequest.GetCuthi(), MyRequest.GetCutinc());
  722.     } catch (NCBI_NS_STD::exception& e) {
  723. ERR_POST(Info << "Exception caught in CSearch::Search: " << e.what());
  724. throw;
  725.     }
  728.     return 0;
  729. }
  732. void CSearch::SetResult(CMSPeakSet& PeakSet, CMSResponse& MyResponse,
  733. double ThreshStart, double ThreshEnd,
  734. double ThreshInc)
  735. {
  736.     CMSPeak* Peaks;
  737.     TPeakSet::iterator iPeakSet;
  739.     TScoreList ScoreList;
  740.     TScoreList::iterator iScoreList;
  741.     CMSHit * MSHit;
  742.     unsigned header;
  743.     char *blastdefline;
  744.     SeqId *sip, *bestid;
  745.     int length;
  746.     unsigned char *Sequence;
  747.     int iSearch; // counter for the length of hit list
  749.     for(iPeakSet = PeakSet.GetPeaks().begin();
  750. iPeakSet != PeakSet.GetPeaks().end();
  751. iPeakSet++ ) {
  752. Peaks = *iPeakSet;
  754. // add to hitset
  755. CRef< CMSHitSet > HitSet (new CMSHitSet);
  756. if(!HitSet) {
  757.     ERR_POST(Error << "omssa: unable to allocate hitset");
  758.     return;
  759. }
  760. HitSet->SetNumber(Peaks->GetNumber());
  761. HitSet->SetFilename(Peaks->GetName());
  762. MyResponse.SetHitsets().push_back(HitSet);
  764. if(Peaks->GetError() == eMSHitError_notenuffpeaks) {
  765.     ERR_POST(Info << "empty set");
  766.     HitSet->SetError(eMSHitError_notenuffpeaks);
  767.     continue;
  768. }
  770. double Threshold, MinThreshold(ThreshStart), MinEval(1000000.0L);
  771. // now calculate scores and sort
  772. for(Threshold = ThreshStart; Threshold <= ThreshEnd; 
  773.     Threshold += ThreshInc) {
  774.     CalcNSort(ScoreList, Threshold, Peaks, false);
  775.     if(!ScoreList.empty()) {
  776. _TRACE("Threshold = " << Threshold <<
  777.        "EVal = " << ScoreList.begin()->first);
  778.     }
  779.     if(!ScoreList.empty() && ScoreList.begin()->first < MinEval) {
  780. MinEval = ScoreList.begin()->first;
  781. MinThreshold = Threshold;
  782.     }
  783.     ScoreList.clear();
  784. }
  785. _TRACE("Min Threshold = " << MinThreshold);
  786. CalcNSort(ScoreList,
  787. #ifdef MSSTATRUN
  788.  ThreshStart
  789. #else
  790.  MinThreshold
  791. #endif
  792. , Peaks, true);
  795. // keep a list of redundant peptides
  796. map <string, CMSHits * > PepDone;
  797. // add to hitset by score
  798. for(iScoreList = ScoreList.begin(), iSearch = 0;
  799.     iScoreList != ScoreList.end();
  800.     iScoreList++, iSearch++) {
  801.     
  802.     CMSHits * Hit;
  803.     CMSPepHit * Pephit;
  805.     MSHit = iScoreList->second;
  806.     header = 0;
  807.     while (readdb_get_header(rdfp, MSHit->GetSeqIndex(), &header,
  808.      &sip,
  809.      &blastdefline)) {
  810. bestid = SeqIdFindBest(sip, SEQID_GI);
  811. if(!bestid) continue;
  813. string seqstring;
  814. string deflinestring(blastdefline);
  815. int iseq;
  816. length = readdb_get_sequence(rdfp, MSHit->GetSeqIndex(),
  817.      &Sequence);
  819. for (iseq = MSHit->GetStart(); iseq <= MSHit->GetStop();
  820.      iseq++) {
  821.     seqstring += UniqueAA[Sequence[iseq]];
  822. }
  823. #ifdef CHECKGI
  824. CheckGi(bestid->data.intvalue);
  825. #endif
  826. if(PepDone.find(seqstring) != PepDone.end()) {
  827.     Hit = PepDone[seqstring];
  828. }
  829. else {
  830.     Hit = new CMSHits;
  831.     Hit->SetPepstring(seqstring);
  832.     double Score = iScoreList->first;
  833.     Hit->SetEvalue(Score);
  834.     Hit->SetPvalue(Score/Peaks->
  835.    GetPeptidesExamined(MSHit->
  836.        GetCharge()));    
  837.     Hit->SetCharge(MSHit->GetCharge());
  838.     CRef<CMSHits> hitref(Hit);
  839.     HitSet->SetHits().push_back(hitref);  
  840.     PepDone[seqstring] = Hit;
  841. }
  843. Pephit = new CMSPepHit;
  845. Pephit->SetStart(MSHit->GetStart());
  846. Pephit->SetStop(MSHit->GetStop());
  847. Pephit->SetGi(bestid->data.intvalue);
  848. Pephit->SetDefline(deflinestring);
  849. CRef<CMSPepHit> pepref(Pephit);
  850. Hit->SetPephits().push_back(pepref);
  851. MemFree(blastdefline);
  852. SeqIdSetFree(sip);
  853.     }
  854. }
  855. ScoreList.clear();
  856.     }
  857. }
  860. // calculate the evalues of the top hits and sort
  861. void CSearch::CalcNSort(TScoreList& ScoreList, double Threshold, CMSPeak* Peaks, bool NewScore)
  862. {
  863.     int iCharges;
  864.     int iHitList;
  865.     int length;
  866.     unsigned char *Sequence;
  868.     for(iCharges = 0; iCharges < Peaks->GetNumCharges(); iCharges++) {
  870. TMSHitList& HitList = Peaks->GetHitList(iCharges);   
  871. for(iHitList = 0; iHitList != Peaks->GetHitListIndex(iCharges);
  872.     iHitList++) {
  873. #ifdef CHECKGI
  874.     unsigned header;
  875.     char *blastdefline;
  876.     SeqId *sip, *bestid;
  878.     header = 0;
  879.     readdb_get_header(rdfp, HitList[iHitList].GetSeqIndex(),
  880.       &header, &sip,&blastdefline);
  881.     bestid = SeqIdFindBest(sip, SEQID_GI);
  882.     if(!bestid) continue;
  883.     CheckGi(bestid->data.intvalue);
  884.     MemFree(blastdefline);
  885.     SeqIdSetFree(sip);
  886. #endif
  888.     // recompute ladder
  889.     length = readdb_get_sequence(rdfp, 
  890.  HitList[iHitList].GetSeqIndex(),
  891.  &Sequence);
  893.     int tempMass = HitList[iHitList].GetMass();
  894.     int tempStart = HitList[iHitList].GetStart();
  895.     int tempStop = HitList[iHitList].GetStop();
  896.     int Charge = HitList[iHitList].GetCharge();
  897.     int Which = Peaks->GetWhich(Charge);
  899.     double a = CalcPoissonMean(tempStart, tempStop, tempMass, Peaks,
  900.        Charge, 
  901.        Threshold);
  902.     if ( a < 0) continue;  // error return;
  903.     if(HitList[iHitList].GetHits() < a) continue;
  904.     double pval;
  905.     if(NewScore) {
  906. int High, Low, NumPeaks, NumLo, NumHi;
  907. Peaks->HighLow(High, Low, NumPeaks, tempMass, Charge, Threshold, NumLo, NumHi);
  908.  
  909. double TopHitProb = ((double)MSNUMTOP)/NumPeaks;
  910. double Normal = CalcNormalTopHit(a, TopHitProb);
  911. pval = CalcPvalueTopHit(a, 
  912. HitList[iHitList].GetHits(Threshold, Peaks->GetMaxI(Which)),
  913. Peaks->GetPeptidesExamined(Charge), Normal, TopHitProb);
  914. #ifdef MSSTATRUN
  915. cout << pval << " " << a << " " << Charge << endl;
  916. // cerr << "pval = " << pval << " Normal = " << Normal << endl;
  917. #endif
  918.     }
  919.     else {
  920. pval =
  921.     CalcPvalue(a, HitList[iHitList].
  922.        GetHits(Threshold, Peaks->GetMaxI(Which)),
  923.        Peaks->GetPeptidesExamined(Charge));
  924.     }
  925.     double eval = pval * Peaks->GetPeptidesExamined(Charge);
  926.     ScoreList.insert(pair<const double, CMSHit *> 
  927.      (eval, &(HitList[iHitList])));
  928. }   
  929.     } 
  930. }
  933. //calculates the poisson times the top hit probability
  934. double CSearch::CalcPoissonTopHit(double Mean, int i, double TopHitProb)
  935. {
  936. double retval;
  937. #ifdef NCBI_OS_MSWIN
  938.     retval =  exp(-Mean) * pow(Mean, i) / exp(LnGamma(i+1));
  939. #else
  940.     retval =  exp(-Mean) * pow(Mean, i) / exp(lgamma(i+1));
  941. #endif
  942. retval *= 1.0L - pow((1.0L-TopHitProb), i);
  943. return retval;
  944. }
  946. double CSearch::CalcPoisson(double Mean, int i)
  947. {
  948. #ifdef NCBI_OS_MSWIN
  949.     return exp(-Mean) * pow(Mean, i) / exp(LnGamma(i+1));
  950. #else
  951.     return exp(-Mean) * pow(Mean, i) / exp(lgamma(i+1));
  952. #endif
  953. }
  955. double CSearch::CalcPvalue(double Mean, int Hits, int n)
  956. {
  957.     if(Hits <= 0) return 1.0L;
  959.     int i;
  960.     double retval(0.0L);
  962.     for(i = 0; i < Hits; i++)
  963. retval += CalcPoisson(Mean, i);
  965.     retval = 1.0L - pow(retval, n);
  966.     if(retval <= 0.0L) retval = numeric_limits<double>::min();
  967.     return retval;
  968. }
  970. #define MSDOUBLELIMIT 0.0000000000000001L
  971. // do full summation of probability distribution
  972. double CSearch::CalcNormalTopHit(double Mean, double TopHitProb)
  973. {
  974.     int i;
  975.     double retval(0.0L), before(-1.0L), increment;
  977.     for(i = 1; i < 1000; i++) {
  978. increment = CalcPoissonTopHit(Mean, i, TopHitProb);
  979. // convergence hack -- on linux (at least) convergence test doesn't work
  980. // for gcc release build
  981. if(increment <= MSDOUBLELIMIT) break;
  982. // if(increment <= numeric_limits<double>::epsilon()) break;
  983. retval += increment;
  984. if(retval == before) break;  // convergence
  985. before = retval;
  986.     }
  987.     return retval;
  988. }
  991. double CSearch::CalcPvalueTopHit(double Mean, int Hits, int n, double Normal, double TopHitProb)
  992. {
  993.     if(Hits <= 0) return 1.0L;
  995.     int i;
  996.     double retval(0.0L), increment;
  998.     for(i = 1; i < Hits; i++) {
  999. increment = CalcPoissonTopHit(Mean, i, TopHitProb);
  1000. if(increment <= MSDOUBLELIMIT) break;
  1001. retval += increment;
  1002.     }
  1004.     retval /= Normal;
  1006.     retval = 1.0L - pow(retval, n);
  1007.     if(retval <= 0.0L) retval = numeric_limits<double>::min();
  1008.     return retval;
  1009. }
  1012. double CSearch::CalcPoissonMean(int Start, int Stop, int Mass, CMSPeak *Peaks,
  1013. int Charge, double Threshold)
  1014. {
  1015.     double retval(1.0);
  1016.     
  1017.     double m; // mass of ladder
  1018.     double o; // min m/z value of peaks
  1019.     double r; // max m/z value of peaks
  1020.     double h; // number of calculated m/z values in one ion series
  1021.     double v1; // number of m/z values above mh/2 and below mh (+1 product ions)
  1022.     double v2; // number of m/z value below mh/2
  1023.     double v;  // number of peaks
  1024.     double t; // mass tolerance width
  1026.     int High, Low, NumPeaks, NumLo, NumHi;
  1027.     Peaks->HighLow(High, Low, NumPeaks, Mass, Charge, Threshold, NumLo, NumHi);
  1028.     if(High == -1) return -1.0; // error
  1029.     o = (double)Low/MSSCALE;
  1030.     r = (double)High/MSSCALE;
  1031.     v1 = NumHi;
  1032.     v2 = NumLo;
  1033.     v = NumPeaks;
  1034.     m = Mass/(double)MSSCALE;
  1035.     h = Stop - Start;
  1036.     t = (Peaks->GetTol())/(double)MSSCALE;
  1037.     // see 12/13/02 notebook, pg. 127
  1039.     retval = 4.0 * t * h * v / m;
  1040.     if(Charge >= kConsiderMult) {
  1041. retval *= (m - 3*o + r)/(r - o);
  1042.     }    
  1043. #if 0
  1044.     // variation that counts +1 and +2 separately
  1045.     // see 8/19/03 notebook, pg. 71
  1046.     retval = 4.0 * t * h / m;
  1047.     if(Charge >= kConsiderMult) {
  1048. // retval *= (m - 3*o + r)/(r - o);
  1049. retval *= (3 * v2 + v1);
  1050.     }
  1051.     else retval *= (v2 + v1);
  1052. #endif
  1053.     return retval;
  1054. }
  1058. CSearch::~CSearch()
  1059. {
  1060.     if(rdfp) readdb_destruct(rdfp);
  1061. }
  1063. /*
  1064. $Log: omssa.cpp,v $
  1065. Revision 1000.4  2004/06/01 18:09:10  gouriano
  1066. PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.20
  1068. Revision 1.20  2004/05/27 20:52:15  lewisg
  1069. better exception checking, use of AutoPtr, command line parsing
  1071. Revision 1.19  2004/05/21 21:41:03  gorelenk
  1072. Added PCH ncbi_pch.hpp
  1074. Revision 1.18  2004/04/06 19:53:20  lewisg
  1075. allow adjustment of precursor charges that allow multiply charged product ions
  1077. Revision 1.17  2004/04/05 20:49:16  lewisg
  1078. fix missed mass bug and reorganize code
  1080. Revision 1.16  2004/03/31 02:00:26  ucko
  1081. +<algorithm> for sort()
  1083. Revision 1.15  2004/03/30 19:36:59  lewisg
  1084. multiple mod code
  1086. Revision 1.14  2004/03/16 20:18:54  gorelenk
  1087. Changed includes of private headers.
  1089. Revision 1.13  2004/03/01 18:24:08  lewisg
  1090. better mod handling
  1092. Revision 1.12  2003/12/22 21:58:00  lewisg
  1093. top hit code and variable mod fixes
  1095. Revision 1.11  2003/12/04 23:39:08  lewisg
  1096. no-overlap hits and various bugfixes
  1098. Revision 1.10  2003/11/20 22:32:50  lewisg
  1099. fix end of sequence X bug
  1101. Revision 1.9  2003/11/20 15:40:53  lewisg
  1102. fix hitlist bug, change defaults
  1104. Revision 1.8  2003/11/17 18:36:56  lewisg
  1105. fix cref use to make sure ref counts are decremented at loop termination
  1107. Revision 1.7  2003/11/14 20:28:06  lewisg
  1108. scale precursor tolerance by charge
  1110. Revision 1.6  2003/11/13 19:07:38  lewisg
  1111. bugs: iterate completely over nr, don't initialize blastdb by iteration
  1113. Revision 1.5  2003/11/10 22:24:12  lewisg
  1114. allow hitlist size to vary
  1116.   Revision 1.4  2003/10/24 21:28:41  lewisg
  1117.   add omssa, xomssa, omssacl to win32 build, including dll
  1118.   
  1119. Revision 1.3  2003/10/22 15:03:32  lewisg
  1120. limits and string compare changed for gcc 2.95 compatibility
  1122.   Revision 1.2  2003/10/21 21:12:17  lewisg
  1123.   reorder headers
  1124.   
  1125. Revision 1.1  2003/10/20 21:32:13  lewisg
  1126. ommsa toolkit version
  1128.   Revision 1.18  2003/10/07 18:02:28  lewisg
  1129.   prep for toolkit
  1130.   
  1131. Revision 1.17  2003/10/06 18:14:17  lewisg
  1132. threshold vary
  1134.   Revision 1.16  2003/08/14 23:49:22  lewisg
  1135.   first pass at variable mod
  1136.   
  1137. Revision 1.15  2003/08/06 18:29:11  lewisg
  1138. support for filenames, numbers using regex
  1140.   Revision 1.14  2003/07/21 20:25:03  lewisg
  1141.   fix missing peak bug
  1142.   
  1143. Revision 1.13  2003/07/19 15:07:38  lewisg
  1144. indexed peaks
  1146.   Revision 1.12  2003/07/18 20:50:34  lewisg
  1147.   *** empty log message ***
  1148.   
  1149. Revision 1.11  2003/07/17 18:45:50  lewisg
  1150. multi dta support
  1152.   Revision 1.10  2003/07/07 16:17:51  lewisg
  1153.   new poisson distribution and turn off histogram
  1154.   
  1155. Revision 1.9  2003/05/01 14:52:10  lewisg
  1156. fixes to scoring
  1158.   Revision 1.8  2003/04/18 20:46:52  lewisg
  1159.   add graphing to omssa
  1160.   
  1161. Revision 1.7  2003/04/07 16:47:16  lewisg
  1162. pc fixes
  1164.   Revision 1.6  2003/04/04 18:43:51  lewisg
  1165.   tax cut
  1166.   
  1167. Revision 1.5  2003/04/02 18:49:51  lewisg
  1168. improved score, architectural changes
  1170.   Revision 1.4  2003/03/21 21:14:40  lewisg
  1171.   merge ming's code, other stuff
  1172.   
  1173. Revision 1.3  2003/02/10 19:37:56  lewisg
  1174. perf and web page cleanup
  1176.   Revision 1.2  2003/02/07 16:18:23  lewisg
  1177.   bugfixes for perf and to work on exp data
  1178.   
  1179. Revision 1.1  2003/02/03 20:39:02  lewisg
  1180. omssa cgi
  1182.   
  1184. */