omssa.cpp

ncbi源码

/*
 * ===========================================================================
 * PRODUCTION $Log: omssa.cpp,v $
 * PRODUCTION Revision 1000.4  2004/06/01 18:09:10  gouriano
 * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.20
 * PRODUCTION
 * ===========================================================================
 */

/* $Id: omssa.cpp,v 1000.4 2004/06/01 18:09:10 gouriano Exp $
* ===========================================================================
*
*                            PUBLIC DOMAIN NOTICE
*               National Center for Biotechnology Information
*
*  This software/database is a "United States Government Work" under the
*  terms of the United States Copyright Act.  It was written as part of
*  the author's official duties as a United States Government employee and
*  thus cannot be copyrighted.  This software/database is freely available
*  to the public for use. The National Library of Medicine and the U.S.
*  Government have not placed any restriction on its use or reproduction.
*
*  Although all reasonable efforts have been taken to ensure the accuracy
*  and reliability of the software and data, the NLM and the U.S.
*  Government do not and cannot warrant the performance or results that
*  may be obtained by using this software or data. The NLM and the U.S.
*  Government disclaim all warranties, express or implied, including
*  warranties of performance, merchantability or fitness for any particular
*  purpose.
*
*  Please cite the authors in any work or product based on this material.
*
* ===========================================================================
*
* Author:  Lewis Y. Geer
*  
* File Description:
*    code to do the ms/ms search and score matches
*
* ===========================================================================
*/

#include <ncbi_pch.hpp>
#include <corelib/ncbistd.hpp>
#include <corelib/ncbiargs.hpp>
#include <corelib/ncbiapp.hpp>
#include <corelib/ncbienv.hpp>
#include <corelib/ncbidiag.hpp>
#include <corelib/ncbi_limits.hpp>

#include <fstream>
#include <string>
#include <list>
#include <deque>
#include <algorithm>
#include <math.h>

#include "SpectrumSet.hpp"
#include "omssa.hpp"

#include "nrutil.h"
#include <ncbimath.h>

USING_NCBI_SCOPE;
USING_SCOPE(objects);
USING_SCOPE(omssa);

/////////////////////////////////////////////////////////////////////////////
//
//  CSearch::
//
//  Performs the ms/ms search
//


CSearch::CSearch(void): rdfp(0)
{
}

int CSearch::InitBlast(const char *blastdb, bool InitDb)
{
    if(!blastdb) return 0;
    if(rdfp) readdb_destruct(rdfp);
    rdfp = readdb_new_ex2((char *)blastdb,TRUE, READDB_NEW_DO_TAXDB | 
			  READDB_NEW_INDEX, NULL, NULL);
    if(!rdfp) return 1;
	
    numseq = readdb_get_num_entries(rdfp);
	
    if(InitDb) {
	int iSearch, length;
	char testchar;
	unsigned char *Sequence;
		
	for(iSearch = 0; iSearch < numseq; iSearch++) {
	    length = readdb_get_sequence(rdfp, iSearch, &Sequence);
	    testchar = Sequence[0];
	}
    }
    return 0;
	
}


// create the ladders from sequence

int CSearch::CreateLadders(unsigned char *Sequence, int iSearch, int position,
			   int endposition,
			   int *Masses, int iMissed, CAA& AA, 
			   CLadder& BLadder,
			   CLadder& YLadder, CLadder& B2Ladder,
			   CLadder& Y2Ladder,
			   unsigned ModMask,
			   const char **Site,
			   int *DeltaMass,
			   int NumMod)
{
    if(!BLadder.CreateLadder(kBIon, 1, (char *)Sequence, iSearch,
			     position, endposition, Masses[iMissed], 
			     MassArray, AA, ModMask, Site, DeltaMass,
			     NumMod)) return 1;
    if(!YLadder.CreateLadder(kYIon, 1, (char *)Sequence, iSearch,
			 position, endposition, Masses[iMissed], 
			     MassArray, AA, ModMask, Site, DeltaMass,
			     NumMod)) return 1;
    B2Ladder.CreateLadder(kBIon, 2, (char *)Sequence, iSearch,
			  position, endposition, 
			  Masses[iMissed], 
			  MassArray, AA, ModMask, Site, DeltaMass,
			  NumMod);
    Y2Ladder.CreateLadder(kYIon, 2, (char *)Sequence, iSearch,
			  position, endposition,
			  Masses[iMissed], 
			  MassArray, AA, ModMask, Site, DeltaMass,
			  NumMod);
    
    return 0;
}


// compare ladders to experiment
int CSearch::CompareLadders(CLadder& BLadder,
			    CLadder& YLadder, CLadder& B2Ladder,
			    CLadder& Y2Ladder, CMSPeak *Peaks,
			    bool OrLadders,  TMassPeak *MassPeak)
{
    if(MassPeak && MassPeak->Charge >= kConsiderMult) {
	Peaks->CompareSorted(BLadder, MSCULLED2, 0); 
	Peaks->CompareSorted(YLadder, MSCULLED2, 0); 
	Peaks->CompareSorted(B2Ladder, MSCULLED2, 0); 
	Peaks->CompareSorted(Y2Ladder, MSCULLED2, 0);
	if(OrLadders) {
	    BLadder.Or(B2Ladder);
	    YLadder.Or(Y2Ladder);	
	}
    }
    else {
	Peaks->CompareSorted(BLadder, MSCULLED1, 0); 
	Peaks->CompareSorted(YLadder, MSCULLED1, 0); 
    }
    return 0;
}

// compare ladders to experiment
bool CSearch::CompareLaddersTop(CLadder& BLadder,
			    CLadder& YLadder, CLadder& B2Ladder,
			    CLadder& Y2Ladder, CMSPeak *Peaks,
			    TMassPeak *MassPeak)
{
    if(MassPeak && MassPeak->Charge >=  kConsiderMult ) {
	if(Peaks->CompareTop(BLadder)) return true; 
	if(Peaks->CompareTop(YLadder)) return true; 
	if(Peaks->CompareTop(B2Ladder)) return true; 
	if(Peaks->CompareTop(Y2Ladder)) return true;
    }
    else {
	if(Peaks->CompareTop(BLadder)) return true; 
	if(Peaks->CompareTop(YLadder)) return true; 
    }
    return false;
}

#ifdef _DEBUG
#define CHECKGI
#ifdef CHECKGI
bool CheckGi(int gi)
{
    if(gi == 21742354) {
	ERR_POST(Info << "test seq");
	return true;
    }
    return false;
}
#endif
#endif

// loads spectra into peaks
void CSearch::Spectrum2Peak(CMSRequest& MyRequest, CMSPeakSet& PeakSet)
{
    CSpectrumSet::Tdata::const_iterator iSpectrum;
    CMSPeak* Peaks;
	
    iSpectrum = MyRequest.GetSpectra().Get().begin();
    for(; iSpectrum != MyRequest.GetSpectra().Get().end(); iSpectrum++) {
	CRef <CMSSpectrum> Spectrum =  *iSpectrum;
	if(!Spectrum) {
	    ERR_POST(Error << "omssa: unable to find spectrum");
	    return;
	}
		
	Peaks = new CMSPeak(MyRequest.GetHitlistlen());
	if(!Peaks) {
	    ERR_POST(Error << "omssa: unable to allocate CMSPeak");
	    return;
	}
	if(Peaks->Read(*Spectrum, MyRequest.GetMsmstol()) != 0) {
	    ERR_POST(Error << "omssa: unable to read spectrum into CMSPeak");
	    return;
	}
	if(Spectrum->CanGetName()) Peaks->SetName() = Spectrum->GetName();
	if(Spectrum->CanGetNumber())
	    Peaks->SetNumber() = Spectrum->GetNumber();
		
	Peaks->Sort();
	Peaks->SetComputedCharge(3);
	Peaks->InitHitList();
	Peaks->CullAll(MyRequest.GetCutlo(), MyRequest.GetSinglewin(),
		       MyRequest.GetDoublewin(), MyRequest.GetSinglenum(),
		       MyRequest.GetDoublenum());
		
	if(Peaks->GetNum(MSCULLED1) < MSPEAKMIN)  {
	    ERR_POST(Info << "omssa: not enough peaks in spectra");
	    Peaks->SetError(eMSHitError_notenuffpeaks);
	}
	PeakSet.AddPeak(Peaks);
		
    }
    PeakSet.SortPeaks(static_cast <int> (MyRequest.GetPeptol()*MSSCALE));
	
}

// compares TMassMasks.  Lower m/z first in sort.
struct CMassMaskCompare {
    bool operator() (const TMassMask& x, const TMassMask& y)
    {
	if(x.Mass < y.Mass) return true;
	return false;
    }
};

// update sites and masses for new peptide
void CSearch::UpdateWithNewPep(int Missed,
			       const char *PepStart[],
			       const char *PepEnd[], 
			       int NumMod[], 
			       const char *Site[][MAXMOD],
			       int DeltaMass[][MAXMOD],
			       int Masses[],
			       int EndMasses[])
{
    // iterate over missed cleavages
    int iMissed;
    // maximum mods allowed
    int ModMax; 
    // iterate over mods
    int iMod;
    

    // update the longer peptides to add the new peptide (Missed-1) on the end
    for(iMissed = 0; iMissed < Missed - 1; iMissed++) {
	// skip start
	if(PepStart[iMissed] == (const char *)-1) continue;
	// reset the end sequences
	PepEnd[iMissed] = PepEnd[Missed - 1];
				
	// update new mod masses to add in any new mods from new peptide

	// first determine the maximum value for updated mod list
	if(NumMod[iMissed] + NumMod[Missed-1] >= MAXMOD)
	    ModMax = MAXMOD - NumMod[iMissed];
	else ModMax = NumMod[Missed-1];

	// now interatu thru the new entries
	for(iMod = NumMod[iMissed]; 
	    iMod < NumMod[iMissed] + ModMax;
	    iMod++) {
	    Site[iMissed][iMod] = 
		Site[Missed-1][iMod - NumMod[iMissed]];

	    DeltaMass[iMissed][iMod] = 
		DeltaMass[Missed-1][iMod - NumMod[iMissed]];

	    //		    Masses[iMissed][iMod] = 
	    //			Masses[Missed-1][iMod - NumMod[iMissed] + 1] +
	    //			Masses[iMissed][NumMod[iMissed] - 1];
	}
				
	// update old masses
	Masses[iMissed] += Masses[Missed - 1];
				
	// update end masses
	EndMasses[iMissed] = EndMasses[Missed - 1];
			
	// update number of Mods
	NumMod[iMissed] += ModMax;
    }	    
}


// create the various combinations of mods
void CSearch::CreateModCombinations(int Missed,
				    const char *PepStart[],
				    TMassMask MassAndMask[][MAXMOD2],
				    int Masses[],
				    int EndMasses[],
				    int NumMod[],
				    int DeltaMass[][MAXMOD],
				    unsigned NumMassAndMask[])
{
    // need to iterate thru combinations that have iMod.
    // i.e. iMod = 3 and NumMod=5
    // 00111, 01011, 10011, 10101, 11001, 11010, 11100, 01101,
    // 01110
    // i[0] = 0 --> 5-3, i[1] = i[0]+1 -> 5-2, i[3] = i[1]+1 -> 5-1
    // then construct bool mask
	    
    // holders for calculated modification mask and modified peptide masses
    unsigned Mask, MassOfMask;
    // iterate thru active mods
    int iiMod;
    // keep track of the number of unique masks created.  each corresponds to a ladder
    int iModCount;
    // missed cleavage
    int iMissed;
    // number of mods to consider
    int iMod;
    // positions of mods
    int ModIndex[MAXMOD];

    // go thru missed cleaves
    for(iMissed = 0; iMissed < Missed; iMissed++) {
	// skip start
	if(PepStart[iMissed] == (const char *)-1) continue;
	iModCount = 0;

	// set up non-modified mass
	MassAndMask[iMissed][iModCount].Mass = 
	    Masses[iMissed] + EndMasses[iMissed];
	MassAndMask[iMissed][iModCount].Mask = 0;
	iModCount++;
		
	// go thru number of mods allowed
	for(iMod = 0; iMod < NumMod[iMissed] && iModCount < MAXMOD2; iMod++) {
		    
	    // initialize ModIndex that points to mod sites
	    InitModIndex(ModIndex, iMod);
	    do {
			
		// calculate mass
		MassOfMask = Masses[iMissed] + EndMasses[iMissed];
		for(iiMod = 0; iiMod <= iMod; iiMod++ ) 
		    MassOfMask += DeltaMass[iMissed][ModIndex[iiMod]];
		// make bool mask
		Mask = MakeBoolMask(ModIndex, iMod);
		// put mass and mask into storage
		MassAndMask[iMissed][iModCount].Mass = MassOfMask;
		MassAndMask[iMissed][iModCount].Mask = Mask;
		// keep track of the  number of ladders
		iModCount++;

	    } while(iModCount < MAXMOD2 &&
		    CalcModIndex(ModIndex, iMod, NumMod[iMissed]));
	} // iMod

	// sort mask and mass by mass
	sort(MassAndMask[iMissed], MassAndMask[iMissed] + iModCount,
	     CMassMaskCompare());
	// keep track of number of MassAndMask
	NumMassAndMask[iMissed] = iModCount;

    } // iMissed
}

//#define MSSTATRUN

int CSearch::Search(CMSRequest& MyRequest, CMSResponse& MyResponse)
{

    try {

	int length;
	CTrypsin trypsin;
	unsigned header;

#ifdef CHECKGI