mspeak.cpp

ncbi源码

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

/* $Id: mspeak.cpp,v 1000.4 2004/06/01 18:09:04 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 deal with spectra and m/z ladders
 *
 * ===========================================================================
 */

#include <ncbi_pch.hpp>
#include <corelib/ncbistd.hpp>
#include <corelib/ncbi_limits.h>

#include <algorithm>

#include "mspeak.hpp"


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

/////////////////////////////////////////////////////////////////////////////
//
//  CMSHit::
//
//  Used by CMSPeak class to hold hits
//


// helper function for RecordHits that scans thru a single ladder
void CMSHit::RecordMatchesScan(CLadder& Ladder, int& iHitInfo, CMSPeak *Peaks,
			       int Which)
{
    try {
	TIntensity Intensity(new unsigned [Ladder.size()]);
	Peaks->CompareSorted(Ladder, Which, &Intensity);

	// examine hits array
	unsigned i;
	for(i = 0; i < Ladder.size(); i++) {
	    // if hit, add to hitlist
	    if(Ladder.GetHit()[i] > 0) {
		SetHitInfo(iHitInfo).SetCharge() = (char) Ladder.GetCharge();
		SetHitInfo(iHitInfo).SetIon() = (char) Ladder.GetType();
		SetHitInfo(iHitInfo).SetNumber() = (short) i;
		SetHitInfo(iHitInfo).SetIntensity() = *(Intensity.get() + i);
		iHitInfo++;
	    }
	}
    } catch (NCBI_NS_STD::exception& e) {
	ERR_POST(Info << "Exception caught in CMSHit::RecordMatchesScan: " << e.what());
	throw;
    }

}

// make a record of the ions hit
void CMSHit::RecordMatches(CLadder& BLadder, CLadder& YLadder,
			   CLadder& B2Ladder,
			   CLadder& Y2Ladder, CMSPeak *Peaks)
{
    // create hitlist.  note that this is deleted in the copy operator
    HitInfo.reset(new CMSHitInfo[Hits]);
    // increment thru hithist
    int iHitInfo(0); 
    int Which = Peaks->GetWhich(Charge);

    // scan thru each ladder
    if(Charge >= kConsiderMult) {
	RecordMatchesScan(BLadder, iHitInfo, Peaks, Which);
	RecordMatchesScan(YLadder, iHitInfo, Peaks, Which);
	RecordMatchesScan(B2Ladder, iHitInfo, Peaks, Which);
	RecordMatchesScan(Y2Ladder, iHitInfo, Peaks, Which);
    }
    else {
	RecordMatchesScan(BLadder, iHitInfo, Peaks, Which);
	RecordMatchesScan(YLadder, iHitInfo, Peaks, Which);
    }

}

// return number of hits above threshold
int CMSHit::GetHits(double Threshold, int MaxI)
{
    int i, retval(0);

    for(i = 0; i < Hits; i++)
        if(SetHitInfo(i).GetIntensity() > MaxI*Threshold)
            retval++;
    return retval;
}



/////////////////////////////////////////////////////////////////////////////
//
//  Comparison functions used in CMSPeak for sorting
//

// compares m/z.  Lower m/z first in sort.
struct CMZICompare {
    bool operator() (CMZI x, CMZI y)
    {
	if(x.MZ < y.MZ) return true;
	return false;
    }
};

// compares m/z.  High m/z first in sort.
struct CMZICompareHigh {
    bool operator() (CMZI x, CMZI y)
    {
	if(x.MZ > y.MZ) return true;
	return false;
    }
};

// compares intensity.  Higher intensities first in sort.
struct CMZICompareIntensity {
    bool operator() (CMZI x, CMZI y)
    {
	if(x.Intensity > y.Intensity) return true;
	return false;
    }
};



/////////////////////////////////////////////////////////////////////////////
//
//  CMSPeak::
//
//  Class used to hold spectra and convert to mass ladders
//


void CMSPeak::xCMSPeak(void)
{
    AAMap = AA.GetMap();
    Sorted[0] = Sorted[1] = false;
    MZI[MSCULLED1] = 0;
    MZI[MSCULLED2] = 0;
    MZI[MSORIGINAL] = 0;
    MZI[MSTOPHITS] = 0;
    Used[MSCULLED1] = 0;
    Used[MSCULLED2] = 0;
    Used[MSORIGINAL] = 0;    
    Used[MSTOPHITS] = 0;    
    PlusOne = 0.8L;
    ComputedCharge = eChargeUnknown;
    Error = eMSHitError_none;
}

CMSPeak::CMSPeak(void)
{ 
    xCMSPeak(); 
}

CMSPeak::CMSPeak(int HitListSizeIn): HitListSize(HitListSizeIn)
{  
    xCMSPeak();
}

CMSPeak::~CMSPeak(void)
{  
    delete [] MZI[MSCULLED1];
    delete [] Used[MSCULLED1];
    delete [] MZI[MSCULLED2];
    delete [] Used[MSCULLED2];
    delete [] MZI[MSORIGINAL];
    delete [] Used[MSTOPHITS];
    delete [] MZI[MSTOPHITS];
    int iCharges;
    for(iCharges = 0; iCharges < GetNumCharges(); iCharges++)
	delete [] HitList[iCharges];
}

// add hit to hitlist.  returns true if successful
bool CMSPeak::AddHit(CMSHit& in, CMSHit *& out)
{
    out = 0;
    // initialize index using charge
    int Index = in.GetCharge() - Charges[0];
    // check to see if hitlist is full
    if(HitListIndex[Index] >= HitListSize) {
	// if less or equal hits than recorded min, don't bother
	if(in.GetHits() <= LastHitNum[Index]) return false;  
	int i, min(HitList[Index][0].GetHits()); //the minimum number of hits in the list
	int minpos(0);     // the position of min
	// find the minimum in the hitslist
	for(i = 1; i < HitListSize; i++) {
	    if(min > HitList[Index][i].GetHits()) {
		min = HitList[Index][i].GetHits();
		minpos = i;
	    }
	}
	// keep record of min
	LastHitNum[Index] = min;	
	// replace in list
	HitList[Index][minpos] = in;
	out =  & (HitList[Index][minpos]);
	return true;
    } 
    else {
	// add to end of list
	HitList[Index][HitListIndex[Index]] = in;
	out = & (HitList[Index][HitListIndex[Index]]);
	HitListIndex[Index]++;
	return true;
    }
}


void CMSPeak::AddTotalMass(int massin, int tolin)
{
    TotalMass = massin;
    tol = tolin;
}


void CMSPeak::Sort(int Which)
{
    if(Which < MSORIGINAL || Which >=  MSNUMDATA) return;
    sort(MZI[Which], MZI[Which] + Num[Which], CMZICompare());
    Sorted[Which] = true;
}

bool CMSPeak::Contains(int value, int Which)
{
    CMZI precursor;
    precursor.MZ = value -  tol; // /2;
    CMZI *p = lower_bound(MZI[Which], MZI[Which] + Num[Which], precursor, CMZICompare());
    if(p == MZI[Which] + Num[Which]) return false;
    if(p->MZ < value + tol/*/2*/) return true;
    return false;
}

		
bool CMSPeak::ContainsFast(int value, int Which)
{

    int x, l, r;
    
    l = 0;
    r = Num[Which] - 1;

    while(l <= r) {
        x = (l + r)/2;
        if (MZI[Which][x].MZ < value - tol) 
	    l = x + 1;
        else if (MZI[Which][x].MZ > value + tol)
	    r = x - 1;
	else return true;
    } 
    
    if (x < Num[Which] - 1 && 
	MZI[Which][x+1].MZ < value + tol && MZI[Which][x+1].MZ > value - tol) 
	return true;
    return false;
}

// compare assuming all lists are sorted
// the intensity array holds the intensity if there is a match to the ladder
// returns total number of matches, which may be more than is recorded in the ladder due to overlap
int CMSPeak::CompareSorted(CLadder& Ladder, int Which, TIntensity* Intensity)
{
    unsigned i(0), j(0);
    int retval(0);
    if(Ladder.size() == 0 ||  Num[Which] == 0) return 0;

    do {
	if (MZI[Which][j].MZ < Ladder[i] - tol) {
	    j++;
	    if(j >= Num[Which]) break;
	    continue;
	}
	else if(MZI[Which][j].MZ > Ladder[i] + tol) {
	    i++;
	    if(i >= Ladder.size()) break;
	    continue;
	}
	else {
	    if(Used[Which][j] == 0) {
		Used[Which][j] = 1;
		Ladder.GetHit()[i] = Ladder.GetHit()[i] + 1;
	    }
	    retval++;
	    // record the intensity if requested, used for auto adjust
	    if(Intensity) {
		*(Intensity->get() + i) = MZI[Which][j].Intensity;
	    }
	    j++;
	    if(j >= Num[Which]) break;
	    i++;
	    if(i >= Ladder.size()) break;
	}
    } while(true);
    return retval;
}

int CMSPeak::Compare(CLadder& Ladder, int Which)
{
    unsigned i;
    int retval(0);
    for(i = 0; i < Ladder.size(); i++) {
	if(ContainsFast(Ladder[i], Which)) {
	    retval++;
	    Ladder.GetHit()[i] = Ladder.GetHit()[i] + 1;
	}
    }
    return retval;
}

// compares the top peaks to the ladder.  returns immediately when finding a hit
bool CMSPeak::CompareTop(CLadder& Ladder)
{
    unsigned i;
    for(i = 0; i < Num[MSTOPHITS]; i++) {
	if(Ladder.ContainsFast(MZI[MSTOPHITS][i].MZ, tol)) return true;
    }
    return false;
}


// read in an asn.1 spectrum and initialize peak values
int CMSPeak::Read(CMSSpectrum& Spectrum, double MSMSTolerance)
{
    try {
	int Scale = Spectrum.GetScale();
	TotalMass = Spectrum.GetPrecursormz()*MSSCALE/Scale;
	SetTolerance(MSMSTolerance);
	Charge = *(Spectrum.GetCharge().begin());
	Num[MSORIGINAL] = 0;   
    
	const CMSSpectrum::TMz& Mz = Spectrum.GetMz();
	const CMSSpectrum::TAbundance& Abundance = Spectrum.GetAbundance();
	MZI[MSORIGINAL] = new CMZI [Mz.size()];
	Num[MSORIGINAL] = Mz.size();

	int i;
	for(i = 0; i < Num[MSORIGINAL]; i++) {
	    MZI[MSORIGINAL][i].MZ = Mz[i]*MSSCALE/Scale;
	    MZI[MSORIGINAL][i].Intensity = Abundance[i]/Scale;
	}
	Sort(MSORIGINAL);
    } catch (NCBI_NS_STD::exception& e) {
	ERR_POST(Info << "Exception in CMSPeak::Read: " << e.what());
	throw;
    }

    return 0;
}

// return the number of peaks in a range
int CMSPeak::CountRange(double StartFraction, double StopFraction)
{
    CMZI Start, Stop;
    int Precursor = static_cast <int> (TotalMass/Charge + tol/2.0);
    Start.MZ = static_cast <int> (StartFraction * Precursor);
    Stop.MZ = static_cast <int> (StopFraction * Precursor);
    CMZI *LoHit = lower_bound(MZI[MSORIGINAL], MZI[MSORIGINAL] +
			      Num[MSORIGINAL], Start, CMZICompare());
    CMZI *HiHit = upper_bound(MZI[MSORIGINAL], MZI[MSORIGINAL] + 
			      Num[MSORIGINAL], Stop, CMZICompare());

    return HiHit - LoHit;
}

// return the percentage of peaks below and at the precursor
int CMSPeak::PercentBelow(void)
{
    CMZI precursor;
    precursor.MZ = static_cast <int> (TotalMass/Charge + tol/2.0);
    CMZI *Hit = upper_bound(MZI[MSORIGINAL], MZI[MSORIGINAL] + Num[MSORIGINAL],
			    precursor, CMZICompare());
    Hit++;  // go above the peak
    if(Hit >= MZI[MSORIGINAL] + Num[MSORIGINAL]) return Num[MSORIGINAL];
    return Hit-MZI[MSORIGINAL];
}

// is the data charge +1?
bool CMSPeak::IsPlus1(double PercentBelowIn)
{
    if(PercentBelowIn/Num[MSORIGINAL] > PlusOne) return true;
    return false;
}

// takes the ratio, low/high, of two ranges in the spectrum
double CMSPeak::RangeRatio(double Start, double Middle, double Stop)
{
    int Lo = CountRange(Start, Middle);
    int Hi = CountRange(Middle, Stop);
    return (double)Lo/Hi;
}

// calculate the charge
void CMSPeak::SetComputedCharge(int MaxCharge)
{
    if(IsPlus1(PercentBelow())) {
	ComputedCharge = eCharge1;
	Charges[0] = 1;
	NumCharges = 1;
    }
    else {
#if 0
	NumCharges = 1;
	if(RangeRatio(0.5, 1.0, 1.5) > 1.25 ) {
	    Charges[0] = 2;
	    ComputedCharge = eCharge2;
	}
	else if(RangeRatio(1.0, 1.5, 2.0) > 1.25) {
	    Charges[0] = 3;
	    ComputedCharge = eCharge3;
	}
	else {  // assume +4
	    Charges[0] = 4;
	    ComputedCharge = eCharge4;
	}
#endif
	ComputedCharge = eChargeNot1; 
	Charges[0] = 2;
	Charges[1] = 3;
	NumCharges = 2;
    }
}

// initializes arrays used to track hits
void CMSPeak::InitHitList(void)
{
    int iCharges;
    for(iCharges = 0; iCharges < GetNumCharges(); iCharges++) {
	LastHitNum[iCharges] = MSHITMIN - 1;
	HitListIndex[iCharges] = 0;
	HitList[iCharges] = new CMSHit[HitListSize];
	PeptidesExamined[iCharges] = 0;
    }
}

void CMSPeak::xWrite(std::ostream& FileOut, CMZI *Temp, int Num)
{
    FileOut << (double)TotalMass/MSSCALE << " " << Charge << endl;

    int i;
    unsigned Intensity;
    for(i = 0; i < Num; i++) {
        Intensity = Temp[i].Intensity;
        if(Intensity == 0.0) Intensity = 1;  // make Mascot happy
        FileOut << (double)Temp[i].MZ/MSSCALE << " " << 
	    Intensity << endl;
    }
}

void CMSPeak::Write(std::ostream& FileOut, CMSPeak::EFileType FileType, int Which)
{
    if(!FileOut || FileType != CMSPeak::eDTA) return;
    xWrite(FileOut, MZI[Which], Num[Which]);
}