C++ Spectrum::getRTime方法代码示例

//#ifdef _MSC_VER 
//#ifndef MSSINGLESCAN_MAIN
//int mssinglescan_main(int argc, char * argv[] ) {
//#else
//int main(int argc, char *argv[]){
//#endif
//#else
int main(int argc, char *argv[]){
//#endif


	//Here are all the variable we are going to need
	MSReader r;
	Spectrum s;
  int j;

  if(argc==1){
    printf("DESCRIPTION: Reads an MS/MS spectrum from any MSToolkit supported file type and outputs to screen in MS2 format.\n\n");
    printf("USAGE: MSSingleScan [scan number] [file]\n");
    exit(0);
  };

	r.readFile(argv[2],s,atoi(argv[1]));
  if(s.getScanNumber()==0) exit(-1);
  printf("S\t%d\t%d\t%.*f\n",s.getScanNumber(),s.getScanNumber(),2,s.getMZ());
	if(s.getRTime()>0) printf("I\tRTime\t%.*f\n",4,s.getRTime());
	for(j=0;j<s.sizeZ();j++){
		printf("Z\t%d\t%.*f\n",s.atZ(j).z,2,s.atZ(j).mz);
	};

	for(j=0;j<s.size();j++){
		printf("%.4f %.4f\n",s.at(j).mz,s.at(j).intensity);
	};

  return 0;

};

bool CPeptideDatabase::readMascot(char* mascotFile, char* dataFile){
  MSReader r;
  MascotParser p;
  Spectrum s;
  MascotLite m;

  int i;
  int percent=0;

  p.clear();
  p.readText(mascotFile);

  if(p.isDistiller()){
    cout << mascotFile << "\t" << p.size() << " peptides from Distiller, no need for RT lookup." << endl;
  } else if(p.hasRTime()){
    cout << mascotFile << "\t" << p.size() << " peptides have RT already, no need for lookup." << endl;
  } else {
    cout << mascotFile << "\tLooking up retention times for " << p.size() << " Mascot peptides..." << endl;
    cout << percent;
    r.readFile(dataFile,s);
  }

  for(i=0;i<p.size();i++){
    if(p.isDistiller()||p.hasRTime()){ 
      m.rTime=p[i].rTime;
    } else {
      if( (int)i*100/p.size() > percent){
        percent=(int)i*100/p.size();
        cout << "\b\b\b" << percent;
      }

      r.readFile(NULL,s,p[i].scanNum);
      if(s.getScanNumber()==0){
        cout << "Scan not found. Please check that correct data file is used. Stopping Mascot parsing." << endl;
        return false;
      }
      m.rTime=s.getRTime();
    }
    m.charge=p[i].charge;
    m.fileID=fileID;
    m.monoMass=p[i].zeroMass;
    m.scanNum=p[i].scanNum;
    strcpy(m.peptide,p[i].sequence_Long);
    strcpy(m.protein,p[i].gene);
    mascot.push_back(m);
  }
  cout << "  Done!" << endl;

  fileID++;

  return true;
}

//First derivative method taken from CSpecAnalyze, returns base peak intensity of the set
void CNoiseReduction::FirstDerivativePeaks(Spectrum& sp, int winSize){
  int i,j;
  float maxIntensity;
  int bestPeak;
  bool bLastPos;
  Spectrum gp;

	int nextBest;
	double FWHM;
	Peak_T centroid;

  bLastPos=false;
  for(i=0;i<sp.size()-winSize;i++){

    if(sp.at(i).intensity<sp.at(i+winSize).intensity) {
      bLastPos=true;
      continue;
    } else {
      if(bLastPos){
				bLastPos=false;
	
				//find max and add peak
				maxIntensity=0;
				for(j=i;j<i+winSize;j++){
				  if (sp.at(j).intensity>maxIntensity){
				    maxIntensity=sp.at(j).intensity;
				    bestPeak = j;
				  }
				}

				//Best estimate of Gaussian centroid
				//Get 2nd highest point of peak
				if(bestPeak==sp.size()-1){
					nextBest=bestPeak-1;
				} else if(sp.at(bestPeak-1).intensity > sp.at(bestPeak+1).intensity){
					nextBest=bestPeak-1;
				} else {
					nextBest=bestPeak+1;
				}

				//Get FWHM
				FWHM = calcFWHM(sp.at(bestPeak).mz);

				//Calc centroid MZ (in three lines for easy reading)
				centroid.mz = pow(FWHM,2)*log(sp.at(bestPeak).intensity/sp.at(nextBest).intensity);
				centroid.mz /= GC*(sp.at(bestPeak).mz-sp.at(nextBest).mz);
				centroid.mz += (sp.at(bestPeak).mz+sp.at(nextBest).mz)/2;

				//Calc centroid intensity
				centroid.intensity=(float)(sp.at(bestPeak).intensity/exp(-pow((sp.at(bestPeak).mz-centroid.mz)/FWHM,2)*GC));

				//some peaks are funny shaped and have bad gaussian fit.
				//if error is more than 10%, keep existing intensity
				if( fabs((sp.at(bestPeak).intensity - centroid.intensity) / centroid.intensity * 100) > 10 ||
            //not a good check for infinity
            centroid.intensity>999999999999.9 ||
            centroid.intensity < 0 ) {
					centroid.intensity=sp.at(bestPeak).intensity;
				}

				//Hack until I put in mass ranges
				if(centroid.mz<0 || centroid.mz>2000) {
					//do nothing if invalid mz
				} else {
					gp.add(centroid);
				}
				i+=winSize-1;
      }

    }
  }
  
  int scanNumber=sp.getScanNumber();
  int scanNumber2=sp.getScanNumber(false);
  float rTime=sp.getRTime();
  sp = gp;
  sp.setRTime(rTime);
  sp.setScanNumber(scanNumber);
  sp.setScanNumber(scanNumber2,true);

}

//.........这里部分代码省略.........

      while(true){
        posRight++;
        if(posRight>=(int)bs.size()) { //buffer is too short on right, add spectra
          r->readFile(lastFile,ts,bs[bs.size()-1].getScanNumber());
          r->readFile(NULL,ts);
          if(ts.getScanNumber()==0) {
            posRight--;
            break;
          }
          bs.push_back(ts);
          //ts.getRawFilter(cFilter2,256);
          if(ts.getMsLevel()==cs.msLevel) {
            index=posRight;
            break;
          }
        } else {
          //bs[posRight].getRawFilter(cFilter2,256);
          if(bs[posRight].getMsLevel()==cs.msLevel) {
            index=posRight;
            break;
          }
        }
      }
    }

    if(index==-1)  continue;
   
    //ts=bs[index];
    v.push_back(index);

    numScans++;
  }

  //cout << "Still Averaging: " << ps.getScanNumber() << endl;
  //cout << " with: ";
  //for(i=0;i<v.size();i++) cout << bs[v[i]].getScanNumber() << " ";
  //cout << endl;

  //cout << numScans << " " << v.size() << endl;

  //Match peaks between pivot scan and neighbors
  for(i=0;i<(int)v.size();i++) vPos.push_back(0);
  for(i=0;i<(int)ps.size();i++){ //iterate all points
    prec = c * ps.at(i).mz * ps.at(i).mz / 2;
    match=1;

    for(k=0;k<(int)v.size();k++){ //iterate all neighbors
      dif=100000.0;
      //cout << "Checking " << bs[v[k]].getScanNumber() << " pos " << vPos[k] << endl;

      for(j=vPos[k];j<bs[v[k]].size();j++){ //check if point is a match
        dt=fabs(ps.at(i).mz-bs[v[k]].at(j).mz);
        if(dt<=dif) {
          if(dt<prec) {
            ps.at(i).intensity+=bs[v[k]].at(j).intensity;
            vPos[k]=j+1;
            match++;
            break;
          }
          dif=dt;
        } else {
          vPos[k]=j-1;
          break;
        }
      }

    }

    //if data point was not visible across enough scans, set it to 0
		if(match<cs.boxcarFilter && match<(int)v.size()) ps.at(i).intensity=0.0;

  }


  //Average points and apply cutoff
  for(i=0;i<ps.size();i++) {
    ps.at(i).intensity/=numScans;
    sp.add(ps.at(i));
    //if(ps.at(i).intensity>=cutoff) sp.add(ps.at(i));
  }

  sp.setScanNumber(ps.getScanNumber());
  sp.setScanNumber(ps.getScanNumber(true),true);
  sp.setRTime(ps.getRTime());
  sp.setRawFilter(cFilter1);

  //clear unused buffer
  if(posLeft>0){
    while(posLeft>0){
      bs.pop_front();
      posLeft--;
      posA--;
    }
  }

  //cout << "Done averaging" << endl;

  return true;
}