mspeak.cpp

ncbi源码

// counts the number of peaks above % of maximum peak

int CMSPeak::AboveThresh(double Threshold, int Which)
{
    CMZI *MZISort = new CMZI[Num[Which]];
    unsigned i;
    for(i = 0; i < Num[Which]; i++) MZISort[i] = MZI[Which][i];

    // sort so higher intensities first
    sort(MZISort, MZISort+Num[Which], CMZICompareIntensity());

    for(i = 1; i < Num[Which] &&
	    (double)MZISort[i].Intensity/MZISort[0].Intensity > Threshold; i++);
    delete [] MZISort;
    return i;
}


// Truncate the at the precursor if plus one and set charge to 1
// if charge is erronously set to 1, set it to 2

void CMSPeak::TruncatePlus1(void)
{
    int PercentBelowVal = PercentBelow();
    if(IsPlus1(PercentBelowVal)) {
        Num[MSORIGINAL] = PercentBelowVal;
        Charge = 1;
    }
    else if(Charge == 1) Charge = 2;
}

// functions used in SmartCull
// take out original peaks below a threshold.  assumes original peaks sorted by intensity
void CMSPeak::CullBaseLine(double Threshold, CMZI *Temp, int& TempLen)
{
    unsigned iMZI;
    for(iMZI = 0; iMZI < Num[MSORIGINAL] && MZI[MSORIGINAL][iMZI].Intensity > Threshold * MZI[MSORIGINAL][0].Intensity; iMZI++);
    copy(&MZI[MSORIGINAL][0], &MZI[MSORIGINAL][iMZI], Temp);
    TempLen = iMZI;
}


// cull precursors
void CMSPeak::CullPrecursor(CMZI *Temp, int& TempLen, double Precursor)
{
    // chop out precursors
    int iTemp(0), iMZI;
    
    for(iMZI = 0; iMZI < TempLen; iMZI++) { 
	if(Temp[iMZI].MZ > Precursor - tol && 
	    Temp[iMZI].MZ < Precursor + tol) continue;
	
	Temp[iTemp] = Temp[iMZI];
	iTemp++;
    }
    TempLen = iTemp;
}

// iterate thru peaks, deleting ones that pass the test
void CMSPeak::CullIterate(CMZI *Temp, int& TempLen, TMZIbool FCN)
{
    if(!FCN) return;
    int iTemp(0), iMZI, jMZI;
    set <int> Deleted;
    
    // scan for isotopes, starting at highest peak
    for(iMZI = 0; iMZI < TempLen - 1; iMZI++) { 
	if(Deleted.count(iMZI) != 0) continue;
	for(jMZI = iMZI + 1; jMZI < TempLen; jMZI++) { 
	    if(Deleted.count(jMZI) != 0) continue;
	    if((*FCN)(Temp[iMZI], Temp[jMZI], tol)) {
		Deleted.insert(jMZI);
		//		Temp[iMZI].Intensity += Temp[jMZI].Intensity; // keep the intensity
	    }
	}
    }
    
    for(iMZI = 0; iMZI < TempLen; iMZI++) {
	if(Deleted.count(iMZI) == 0) {
	    Temp[iTemp] = Temp[iMZI];
	    iTemp++;
	}
    }
    TempLen = iTemp;
}

// object for looking for isotopes.  true if isotope
static bool FCompareIsotope(const CMZI& x, const CMZI& y, int tol)
{
    if(y.MZ < x.MZ + 2*MSSCALE + tol && y.MZ > x.MZ - tol) return true;
    return false;
}

// cull isotopes using the Markey Method

void CMSPeak::CullIsotope(CMZI *Temp, int& TempLen)
{
    CullIterate(Temp, TempLen, &FCompareIsotope);
}


// function for looking for H2O or NH3.  true if is.
static bool FCompareH2ONH3(const CMZI& x, const CMZI& y, int tol)
{
    if((y.MZ > x.MZ - 18*MSSCALE - tol && y.MZ < x.MZ - 18*MSSCALE + tol ) ||
       (y.MZ > x.MZ - 17*MSSCALE - tol && y.MZ < x.MZ - 17*MSSCALE + tol))
	return true;
    return false;
}

// cull peaks that are water or ammonia loss
// note that this only culls the water or ammonia loss if these peaks have a lesser
// less intensity
void CMSPeak::CullH20NH3(CMZI *Temp, int& TempLen)
{
    // need to start at top of m/z ladder, work way down.  sort
    CullIterate(Temp, TempLen, &FCompareH2ONH3);
}


// use smartcull on all charges
void CMSPeak::CullAll(double Threshold, int SingleWindow,
		      int DoubleWindow, int SingleHit, int DoubleHit)
{    
    int iCharges;
    for(iCharges = 0; iCharges < GetNumCharges(); iCharges++)
	SmartCull(Threshold, GetCharges()[iCharges], SingleWindow,
		  DoubleWindow, SingleHit, DoubleHit);

    // make the high intensity list
    iCharges = GetNumCharges() - 1;
    int Which = GetWhich(GetCharges()[iCharges]);
    CMZI *Temp = new CMZI [Num[Which]]; // temporary holder
    copy(MZI[Which], MZI[Which] + Num[Which], Temp);
    sort(Temp, Temp + Num[Which], CMZICompareIntensity());

    if(Num[Which] > MSNUMTOP)  Num[MSTOPHITS] = MSNUMTOP;
    else  Num[MSTOPHITS] = Num[Which];

    MZI[MSTOPHITS] = new CMZI [Num[MSTOPHITS]]; // holds highest hits
    Used[MSTOPHITS] = new char [Num[MSTOPHITS]];
    ClearUsed(MSTOPHITS);
    Sorted[MSTOPHITS] = false;
    copy(Temp, Temp + Num[MSTOPHITS], MZI[MSTOPHITS]);
    Sort(MSTOPHITS);

    delete [] Temp;
}

// check to see if TestMZ is Diff away from BigMZ
bool CMSPeak::IsAtMZ(int BigMZ, int TestMZ, int Diff, int tol)
{
    if(TestMZ < BigMZ - Diff*MSSCALE + tol && 
       TestMZ > BigMZ - Diff*MSSCALE - tol)
	return true;
    return false;
}

// see if TestMZ can be associated with BigMZ, e.g. water loss, etc.
bool CMSPeak::IsMajorPeak(int BigMZ, int TestMZ, int tol)
{
    if (IsAtMZ(BigMZ, TestMZ, 1, tol) || 
	IsAtMZ(BigMZ, TestMZ, 16, tol) ||
	IsAtMZ(BigMZ, TestMZ, 17, tol) ||
	IsAtMZ(BigMZ, TestMZ, 18, tol)) return true;
    return false;
}

// recursively culls the peaks
void CMSPeak::SmartCull(double Threshold, int Charge, int SingleWindow,
			int DoubleWindow, int SingleNum, int DoubleNum)
{
    //       sort(MZI[MSORIGINAL], MZI[MSORIGINAL] + Num[MSORIGINAL], CMZICompare());
    sort(MZI[MSORIGINAL], MZI[MSORIGINAL] + Num[MSORIGINAL], CMZICompareIntensity());
    Sorted[MSORIGINAL] = false;
    double Precursor;
    int Which = GetWhich(Charge);
    Precursor = GetMass()/(double)Charge;

    int iMZI = 0;  // starting point

    int TempLen(0);
    CMZI *Temp = new CMZI [Num[MSORIGINAL]]; // temporary holder

    // prep the data
    CullBaseLine(Threshold, Temp, TempLen);
#ifdef DEBUG_PEAKS1
    {
	sort(Temp, Temp+TempLen , CMZICompare());
	ofstream FileOut("afterbase.dta");
	xWrite(FileOut, Temp, TempLen);
	sort(Temp, Temp+TempLen , CMZICompareIntensity());
    }
#endif
    CullPrecursor(Temp, TempLen, Precursor);
#ifdef DEBUG_PEAKS1
    {
	sort(Temp, Temp+TempLen , CMZICompare());
	ofstream FileOut("afterprecurse.dta");
	xWrite(FileOut, Temp, TempLen);
	sort(Temp, Temp+TempLen , CMZICompareIntensity());
    }
#endif
    CullIsotope(Temp, TempLen);
#ifdef DEBUG_PEAKS1
    {
	sort(Temp, Temp+TempLen , CMZICompare());
	ofstream FileOut("afteriso.dta");
	xWrite(FileOut, Temp, TempLen);
	sort(Temp, Temp+TempLen , CMZICompareIntensity());
    }
#endif
    // h20 and nh3 loss ignored until this cut can be studied
#if 0
    sort(Temp, Temp + TempLen, CMZICompareHigh());
    CullH20NH3(Temp, TempLen);
    sort(Temp, Temp + TempLen, CMZICompareIntensity());
#endif

    // now do the recursion
    int iTemp(0), jMZI;
    set <int> Deleted;
    map <int, int> MajorPeak; // maps a peak to it's "major peak"
    int HitsAllowed; // the number of hits allowed in a window
    int Window; // the m/z range used to limit the number of peaks
    int HitCount;  // number of nearby peaks
   
    // scan for isotopes, starting at highest peak
    for(iMZI = 0; iMZI < TempLen - 1; iMZI++) { 
	if(Deleted.count(iMZI) != 0) continue;
	HitCount = 0;
	if(Charge <  kConsiderMult || Temp[iMZI].MZ > Precursor) {
	    // if charge 1 region, allow fewer peaks
	    Window = SingleWindow; //27;
	    HitsAllowed = SingleNum;
	}
	else {
	    // otherwise allow a greater number
	    Window = DoubleWindow; //14;
	    HitsAllowed = DoubleNum;
	}
	// go over smaller peaks
	set <int> Considered;
	for(jMZI = iMZI + 1; jMZI < TempLen; jMZI++) { 
	    if(Deleted.count(jMZI) != 0) continue;
	    if(Temp[jMZI].MZ < Temp[iMZI].MZ + Window*MSSCALE + tol &&
	       Temp[jMZI].MZ > Temp[iMZI].MZ - Window*MSSCALE - tol) {
		if(IsMajorPeak(Temp[iMZI].MZ, Temp[jMZI].MZ, tol)) {
		    // link little peak to big peak
		    MajorPeak[jMZI] = iMZI;
		    // ignore for deletion
		    continue;
		}
		// if linked to a major peak, skip
		if(MajorPeak.find(jMZI) != MajorPeak.end()) {
		    if(Considered.find(jMZI) != Considered.end())
			continue;
		    Considered.insert(jMZI);
		    HitCount++;
		    if(HitCount <= HitsAllowed)  continue;
		}
		// if the first peak, skip
		HitCount++;
		if(HitCount <= HitsAllowed)  continue;
		Deleted.insert(jMZI);
	    }
	}
    }
    
    // delete the culled peaks
    for(iMZI = 0; iMZI < TempLen; iMZI++) {
	if(Deleted.count(iMZI) == 0) {
	    Temp[iTemp] = Temp[iMZI];
	    iTemp++;
	}
    }
    TempLen = iTemp;


    // make the array of culled peaks
    if(MZI[Which]) delete [] MZI[Which];
    if(Used[Which]) delete [] Used[Which];
    Num[Which] = TempLen;
    MZI[Which] = new CMZI [TempLen];
    Used[Which] = new char [TempLen];
    ClearUsed(Which);
    copy(Temp, Temp+TempLen, MZI[Which]);
    Sort(Which);
#ifdef DEBUG_PEAKS1
    {
	ofstream FileOut("aftercull.dta");
	xWrite(FileOut, MZI[Which], Num[Which]);
    }
#endif

    delete [] Temp;
}

// return the lowest culled peak and the highest culled peak less than the
// +1 precursor mass
void CMSPeak::HighLow(int& High, int& Low, int& NumPeaks, int PrecursorMass,
		      int Charge, double Threshold, int& NumLo, int& NumHi)
{
    int Which = GetWhich(Charge);

    if(!Sorted[Which]) Sort(Which);
    if(Num[Which] < 2) {
	High = Low = -1;
	NumPeaks = NumLo = NumHi = 0;
	return;
    }

    Low = PrecursorMass;
    High = 0;
    NumPeaks = 0;
    NumLo = 0;
    NumHi = 0;

    int MaxI = GetMaxI(Which);

    unsigned iMZI;
    for(iMZI = 0; iMZI < Num[Which]; iMZI++) {
	if(MZI[Which][iMZI].Intensity > Threshold*MaxI &&
	   MZI[Which][iMZI].MZ <= PrecursorMass) {
	    if(MZI[Which][iMZI].MZ > High) {
		High = MZI[Which][iMZI].MZ;
	    }
	    if(MZI[Which][iMZI].MZ < Low) {
		Low = MZI[Which][iMZI].MZ;
	    }
	    NumPeaks++;
	    if(MZI[Which][iMZI].MZ < PrecursorMass/2.0) NumLo++;
	    else NumHi++;
	}
    }
}

// count the number of putative AA intervals.
// Nodup true means don't count the peaks twice
int CMSPeak::CountAAIntervals(CMassArray& MassArray, bool Nodup, int Which)
{	
    unsigned ipeaks, low;
    int i;
    int PeakCount(0);
    if(!Sorted[Which]) return -1;
    const int *IntMassArray = MassArray.GetIntMass();

    //	list <int> intensity;
    //	int maxintensity = 0;

    for(ipeaks = 0 ; ipeaks < Num[Which] - 1; ipeaks++) {
	//		if(maxintensity < MZI[ipeaks].Intensity) maxintensity = MZI[ipeaks].Intensity;
	low = ipeaks;
	low++;
	//		if(maxintensity < MZI[low].Intensity) maxintensity = MZI[low].Intensity;
	for(; low < Num[Which]; low++) {
	    for(i = 0; i < kNumUniqueAA; i++) {
		if(IntMassArray[i] == 0) continue;  // skip gaps, etc.
		if(MZI[Which][low].MZ- MZI[Which][ipeaks].MZ < IntMassArray[i] + tol/2.0 &&
		   MZI[Which][low].MZ - MZI[Which][ipeaks].MZ > IntMassArray[i] - tol/2.0 ) {	  
		    PeakCount++;
		    //					intensity.push_back(MZI[ipeaks].Intensity);
		    if(Nodup) goto newpeak;
		    else goto oldpeak;
		}
	    }
	oldpeak:
	    i = i;
	}
    newpeak:
	i = i;
    }

    return PeakCount;
}


// return maximum intensity value
int CMSPeak::GetMaxI(int Which)
{
    unsigned Intensity(0), i;
    for(i = 0; i < Num[Which]; i++) {
        if(Intensity < MZI[Which][i].Intensity)
	    Intensity = MZI[Which][i].Intensity;
    }
    return Intensity;
}



/////////////////////////////////////////////////////////////////////////////
//
//  CMSPeakSet::
//
//  Class used to hold sets of CMSPeak and access them quickly
//


CMSPeakSet::~CMSPeakSet() 
{
    while(!PeakSet.empty()) {
	delete *PeakSet.begin();
	PeakSet.pop_front();
    }
}

// compares m/z.  Lower m/z first in sort.
struct CMassPeakCompareHi {
    bool operator() (TMassPeak x, TMassPeak y)
    {
	if(x.Mass + x.Peptol < y.Mass + y.Peptol) return true;
	return false;
    }
};

void CMSPeakSet::SortPeaks(int Peptol)
{
    int iCharges;
    CMSPeak* Peaks;
    TPeakSet::iterator iPeakSet;
    int CalcMass; // the calculated mass
    TMassPeak temp;
    int ptol; // charge corrected mass tolerance

    // first sort
    for(iPeakSet = GetPeaks().begin();
	iPeakSet != GetPeaks().end();
	iPeakSet++ ) {
	Peaks = *iPeakSet;
	// skip empty spectra
	if(Peaks->GetError() == eMSHitError_notenuffpeaks) continue;

	// loop thru possible charges
	for(iCharges = 0; iCharges < Peaks->GetNumCharges(); iCharges++) {
	    // correction for incorrect charge determination.
	    // see 12/13/02 notebook, pg. 135
	    ptol = Peaks->GetCharges()[iCharges] * Peptol;
	    CalcMass = static_cast <int> ((Peaks->GetMass() +
					   Peaks->GetCharge()*kProton*MSSCALE) * 
					  Peaks->GetCharges()[iCharges]/(double)(Peaks->GetCharge()) - 
					  Peaks->GetCharges()[iCharges]*kProton*MSSCALE);
	    temp.Mass = CalcMass;
	    temp.Peptol = ptol;
	    temp.Charge = Peaks->GetCharges()[iCharges];
	    temp.Peak = Peaks;
	    // order by upper bound
	    MassMap.insert(pair <const int, TMassPeak>(temp.Mass + temp.Peptol, temp)); 
	}
    } 

    // then create static array

    ArraySize = MassMap.size();
    MassPeak.reset(new TMassPeak[ArraySize]);

    TMassPeakMap::iterator iMassMap;
    int i(0);
    for(iMassMap = MassMap.begin(); iMassMap != MassMap.end(); iMassMap++, i++) {
	GetMassPeak(i).Mass =
 iMassMap->second.Mass;
	GetMassPeak(i).Peptol = iMassMap->second.Peptol;
	GetMassPeak(i).Peak = iMassMap->second.Peak;
	GetMassPeak(i).Charge = iMassMap->second.Charge;
    }

    MassMap.clear();

}


// Get the first index into the sorted array where the mass
// + tolerance is >= the given calculated mass. 
TMassPeak *CMSPeakSet::GetIndexLo(int Mass)
{
    TMassPeak temp;
    TMassPeak *retval;
    temp.Mass = Mass;
    temp.Peptol = 0;
    // look for first spectrum whose upper mass value + tolerance exceeds calculated mass
    retval = lower_bound(MassPeak.get(), MassPeak.get() + ArraySize, temp, 
			 CMassPeakCompareHi());
    return retval;
}

// get peak for sorted list by index into list
CMSPeak *CMSPeakSet::GetPeak(int Index)
{
    if(Index < 0 || Index >= ArraySize) return 0;
    return GetMassPeak(Index).Peak;
}