C++ SpectrumPtr::cvParam方法代码示例

void printSpectrumList(const SpectrumList& sl, ostream& os)
{
    os << "size: " << sl.size() << endl;

    for (size_t i=0, end=sl.size(); i<end; i++)
    {
        SpectrumPtr spectrum = sl.spectrum(i, false);
        os << spectrum->index << " " 
           << spectrum->id << " "
           << "ms" << spectrum->cvParam(MS_ms_level).value << " "
           << "scanEvent:" << spectrum->scanList.scans[0].cvParam(MS_preset_scan_configuration).value << " "
           << "scanTime:" << spectrum->scanList.scans[0].cvParam(MS_scan_start_time).timeInSeconds() << " "
           << endl;
    }
}

void cat(const char* filename)
{
    // open the data file

    ExtendedReaderList readers; // for vendor Reader support
    MSDataFile msd(filename, &readers);

    // verify that we have a SpectrumList

    if (!msd.run.spectrumListPtr.get())
        throw runtime_error("[mscat] No spectra found.");

    SpectrumList& spectrumList = *msd.run.spectrumListPtr;

    // write header

    const size_t columnWidth = 14;

    cout << "#"
         << setw(columnWidth) << "scanNumber"
         << setw(columnWidth) << "msLevel"
         << setw(columnWidth) << "m/z"
         << setw(columnWidth) << "intensity" << endl;

    // iterate through the spectra in the SpectrumList
    for (size_t i=0, size=spectrumList.size(); i!=size; i++)
    {
        // retrieve the spectrum, with binary data
        const bool getBinaryData = true;
        SpectrumPtr spectrum = spectrumList.spectrum(i, getBinaryData);

        // fill in MZIntensityPair vector for convenient access to binary data
        vector<MZIntensityPair> pairs;
        spectrum->getMZIntensityPairs(pairs);

        // iterate through the m/z-intensity pairs
        for (vector<MZIntensityPair>::const_iterator it=pairs.begin(), end=pairs.end(); it!=end; ++it)
        {
            cout << " "
                 << setw(columnWidth) << spectrum->nativeID
                 << setw(columnWidth) << "ms" + spectrum->cvParam(MS_ms_level).value
                 << setw(columnWidth) << fixed << setprecision(4) << it->mz
                 << setw(columnWidth) << fixed << setprecision(2) << it->intensity << endl;
        }
    }
}

    void writeSpectrumBinary(SpectrumPtr s, int version, bool compress, ostream& os)
    {
      // todo
      // MSnScanInfo header(s);
      // header.writeSpectrumHeader(version, compress, os);

        int scanNum = getScanNumber(s);
        os.write(reinterpret_cast<char *>(&scanNum), sizeIntMSn);
        os.write(reinterpret_cast<char *>(&scanNum), sizeIntMSn); // Yes, there are two

        Precursor& precur = s->precursors[0];
        SelectedIon& si = precur.selectedIons[0];
        double mz = si.cvParam(MS_selected_ion_m_z).valueAs<double>();
        os.write(reinterpret_cast<char *>(&mz), sizeDoubleMSn);

        float rt = (float) s->scanList.scans[0].cvParam(MS_scan_start_time).timeInSeconds();
        os.write(reinterpret_cast<char *>(&rt), sizeFloatMSn);
        
        if (version >= 2)
        {
            float basePeakIntensity = s->cvParam(MS_base_peak_intensity).valueAs<float>();
            os.write(reinterpret_cast<char *>(&basePeakIntensity), sizeFloatMSn);

            double basePeakMZ = s->cvParam(MS_base_peak_m_z).valueAs<double>();
            os.write(reinterpret_cast<char *>(&basePeakMZ), sizeDoubleMSn);
            
            // We don't have this information, but we need to write something,
            // so pad with 0's. (version 2 specific data)
            double conversionFactorA = (double)0;
            os.write(reinterpret_cast<char *>(&conversionFactorA), sizeDoubleMSn); 
            double conversionFactorB = (double)0;
            os.write(reinterpret_cast<char *>(&conversionFactorB), sizeDoubleMSn); 

            double tic = s->cvParam(MS_total_ion_current).valueAs<double>();
            os.write(reinterpret_cast<char *>(&tic), sizeDoubleMSn);

            // TODO
            float ionInjectionTime = (float)0;
            os.write(reinterpret_cast<char *>(&ionInjectionTime), sizeFloatMSn);
        }
        
        vector<int> charges;
        vector<double> masses;
        int numChargeStates = 0;
        BOOST_FOREACH(const SelectedIon& curIon, precur.selectedIons)
        {
            numChargeStates += getChargeStates(curIon, charges, masses);
        }
        os.write(reinterpret_cast<char *>(&numChargeStates), sizeIntMSn);
        
        bool hasAccurateMass = false;
        if (version == 3)
        {
          int numEzStates = 0;
          CVParam massParam = si.cvParam(MS_accurate_mass);
          if (!massParam.empty())
          {
            numEzStates = numChargeStates;
            hasAccurateMass = true;
          }
          os.write(reinterpret_cast<char *>(&numEzStates), sizeIntMSn);
        }

        int numPeaks = (int) s->defaultArrayLength;
        os.write(reinterpret_cast<char *>(&numPeaks), sizeIntMSn);

        // end spectrum header info

        // Write out each charge state and corresponding mass
        for(int i = 0; i < numChargeStates; i++)
        {
            os.write(reinterpret_cast<char *>(&(charges[i])), sizeIntMSn);
            os.write(reinterpret_cast<char *>(&(masses[i])), sizeDoubleMSn);
        }
    
        // if there are accurate masses, write out EZ entries
        if( hasAccurateMass ){
          float blank = 0;  // we don't have rTime or area, pad with zeros
          for(int i=0; i < numChargeStates; i++){
            os.write(reinterpret_cast<char *>(&charges[i]), sizeIntMSn);
            os.write(reinterpret_cast<char *>(&masses[i]), sizeDoubleMSn);
            os.write(reinterpret_cast<char *>(&blank), sizeFloatMSn);
            os.write(reinterpret_cast<char *>(&blank), sizeFloatMSn);
          }
        }

        // Do we need to write compressed m/z, intensity arrays?
        if (compress)
        {
            writeCompressedPeaks(s, os);
        }
        else
        {
            // No need to compress, just write out the arrays
            const BinaryDataArray& mzArray = *s->getMZArray();
            const BinaryDataArray& intensityArray = *s->getIntensityArray();
            for(int i = 0; i < numPeaks; i++)
            {
                double mzPeak = mzArray.data[i];
                os.write(reinterpret_cast<char *>(&mzPeak), sizeDoubleMSn);
//.........这里部分代码省略.........

void RAMPAdapter::Impl::getScanHeader(size_t index, ScanHeaderStruct& result, bool reservePeaks /*= true*/) const
{
    // use previous spectrum if possible
    if (!lastSpectrum.get() || lastSpectrum->index != index)
        lastSpectrum = msd_.run.spectrumListPtr->spectrum(index, reservePeaks);

    SpectrumPtr spectrum = lastSpectrum;

    Scan dummy;
    Scan& scan = spectrum->scanList.scans.empty() ? dummy : spectrum->scanList.scans[0];

    result.seqNum = static_cast<int>(index + 1);
    result.acquisitionNum = getScanNumber(index);
    result.msLevel = spectrum->cvParam(MS_ms_level).valueAs<int>();
    result.peaksCount = static_cast<int>(spectrum->defaultArrayLength);
    result.totIonCurrent = spectrum->cvParam(MS_total_ion_current).valueAs<double>();
    result.retentionTime = scan.cvParam(MS_scan_start_time).timeInSeconds();
    result.basePeakMZ = spectrum->cvParam(MS_base_peak_m_z).valueAs<double>();    
    result.basePeakIntensity = spectrum->cvParam(MS_base_peak_intensity).valueAs<double>();    
    result.collisionEnergy = 0;
    result.ionisationEnergy = spectrum->cvParam(MS_ionization_energy).valueAs<double>();
    result.lowMZ = spectrum->cvParam(MS_lowest_observed_m_z).valueAs<double>();        
    result.highMZ = spectrum->cvParam(MS_highest_observed_m_z).valueAs<double>();        
    result.precursorScanNum = 0;
    result.precursorMZ = 0;
    result.precursorCharge = 0;
    result.precursorIntensity = 0;
    result.compensationVoltage = 0;
    

    std::string filterLine = scan.cvParam(MS_filter_string).value;
    
    size_t found = filterLine.find("cv=");
    
    if (found!=string::npos) {
      filterLine = filterLine.substr(found+3);
      found = filterLine.find_first_of(" ");
      if (found!=string::npos) {
          filterLine = filterLine.substr(0, found);
          result.compensationVoltage = atof(filterLine.c_str());
      }
    }

    result.filterLine = filterLine;



    if (!spectrum->precursors.empty())
    {
        const Precursor& precursor = spectrum->precursors[0];
        result.collisionEnergy = precursor.activation.cvParam(MS_collision_energy).valueAs<double>();
        size_t precursorIndex = msd_.run.spectrumListPtr->find(precursor.spectrumID);

        if (precursorIndex < msd_.run.spectrumListPtr->size())
        {
            const SpectrumIdentity& precursorSpectrum = msd_.run.spectrumListPtr->spectrumIdentity(precursorIndex);
            string precursorScanNumber = id::translateNativeIDToScanNumber(nativeIdFormat_, precursorSpectrum.id);
            
            if (precursorScanNumber.empty()) // unsupported nativeID type
            {
                // assume scanNumber is a 1-based index, consistent with this->index() method
                result.precursorScanNum = precursorIndex+1;
            } 
            else 
            {
                result.precursorScanNum = lexical_cast<int>(precursorScanNumber);
            }
        }
        if (!precursor.selectedIons.empty())
        {
            result.precursorMZ = precursor.selectedIons[0].cvParam(MS_selected_ion_m_z).valueAs<double>();
            if (!result.precursorMZ)
            { // mzML 1.0?
                result.precursorMZ = precursor.selectedIons[0].cvParam(MS_m_z).valueAs<double>();
            }
            result.precursorCharge = precursor.selectedIons[0].cvParam(MS_charge_state).valueAs<int>();
            result.precursorIntensity = precursor.selectedIons[0].cvParam(MS_peak_intensity).valueAs<double>();
        }
    }

    BOOST_STATIC_ASSERT(SCANTYPE_LENGTH > 4);
    memset(result.scanType, 0, SCANTYPE_LENGTH);
    strcpy(result.scanType, "Full"); // default
    if (spectrum->hasCVParam(MS_zoom_scan))
        strcpy(result.scanType, "Zoom");

    result.mergedScan = 0; // TODO 
    result.mergedResultScanNum = 0; // TODO 
    result.mergedResultStartScanNum = 0; // TODO 
    result.mergedResultEndScanNum = 0; // TODO 
    result.filePosition = spectrum->sourceFilePosition; 
}

PWIZ_API_DECL SpectrumPtr SpectrumList_PeakPicker::spectrum(size_t index, bool getBinaryData) const
{
    SpectrumPtr s;
    
    switch (mode_)
    {
        case 1:
            s = dynamic_cast<detail::SpectrumList_Thermo*>(&*inner_)->spectrum(index, getBinaryData, msLevelsToPeakPick_);
            break;

        case 2:
            s = dynamic_cast<detail::SpectrumList_Bruker*>(&*inner_)->spectrum(index, getBinaryData, msLevelsToPeakPick_);
            break;

        case 3:
            s = dynamic_cast<detail::SpectrumList_ABI*>(&*inner_)->spectrum(index, getBinaryData, msLevelsToPeakPick_);
            break;

        case 4:
            s = dynamic_cast<detail::SpectrumList_Agilent*>(&*inner_)->spectrum(index, getBinaryData, msLevelsToPeakPick_);
            break;

        case 5:
            s = dynamic_cast<detail::SpectrumList_ABI_T2D*>(&*inner_)->spectrum(index, getBinaryData, msLevelsToPeakPick_);
            break;

        case 0:
        default:
            s = inner_->spectrum(index, true);
            break;
    }

    if (!msLevelsToPeakPick_.contains(s->cvParam(MS_ms_level).valueAs<int>()))
        return s;

    vector<CVParam>& cvParams = s->cvParams;
    vector<CVParam>::iterator itr = std::find(cvParams.begin(), cvParams.end(), MS_profile_spectrum);

    // return non-profile spectra as-is
    // (could have been acquired as centroid, or vendor may have done the centroiding)
    if (itr == cvParams.end())
        return s;

    // make sure the spectrum has binary data
    if (!s->getMZArray().get() || !s->getIntensityArray().get())
        s = inner_->spectrum(index, true);

    // replace profile term with centroid term
    vector<CVParam>& cvParams2 = s->cvParams;
    *(std::find(cvParams2.begin(), cvParams2.end(), MS_profile_spectrum)) = MS_centroid_spectrum;

    try
    {
        vector<double>& mzs = s->getMZArray()->data;
        vector<double>& intensities = s->getIntensityArray()->data;
        vector<double> xPeakValues, yPeakValues;
        algorithm_->detect(mzs, intensities, xPeakValues, yPeakValues);
        mzs.swap(xPeakValues);
        intensities.swap(yPeakValues);
        s->defaultArrayLength = mzs.size();
    }
    catch(std::exception& e)
    {
        throw std::runtime_error(std::string("[SpectrumList_PeakPicker::spectrum()] Error picking peaks: ") + e.what());
    }

    s->dataProcessingPtr = dp_;
    return s;
}

PWIZ_API_DECL void ThresholdFilter::operator () (const SpectrumPtr s) const
{
    if (!msLevelsToThreshold.contains(s->cvParam(MS_ms_level).valueAs<int>()))
        return;

    // do nothing to empty spectra
    if (s->defaultArrayLength == 0)
        return;

    if (byType == ThresholdingBy_Count ||
        byType == ThresholdingBy_CountAfterTies)
    {
        // if count threshold is greater than number of data points, return as is
        if (s->defaultArrayLength <= threshold)
            return;
        else if (threshold == 0)
        {
            s->getMZArray()->data.clear();
            s->getIntensityArray()->data.clear();
            s->defaultArrayLength = 0;
            return;
        }
    }

    vector<MZIntensityPair> mzIntensityPairs;
    s->getMZIntensityPairs(mzIntensityPairs);

    if (orientation == Orientation_MostIntense)
        sort(mzIntensityPairs.begin(), mzIntensityPairs.end(), orientationMore_Predicate);
    else if (orientation == Orientation_LeastIntense)
        sort(mzIntensityPairs.begin(), mzIntensityPairs.end(), orientationLess_Predicate);
    else
        throw runtime_error("[threshold()] invalid orientation type");

    double tic = accumulate(mzIntensityPairs.begin(), mzIntensityPairs.end(), 0.0, MZIntensityPairIntensitySum());

    if (tic == 0)
    {
        s->getMZArray()->data.clear();
        s->getIntensityArray()->data.clear();
        s->defaultArrayLength = 0;
        return;
    }

    double bpi = orientation == Orientation_MostIntense ? mzIntensityPairs.front().intensity
                                                        : mzIntensityPairs.back().intensity;

    // after the threshold is applied, thresholdItr should be set to the first data point to erase
    vector<MZIntensityPair>::iterator thresholdItr;

    switch (byType)
    {
        case ThresholdingBy_Count:
            // no need to check bounds on thresholdItr because it gets checked above
            thresholdItr = mzIntensityPairs.begin() + (size_t) threshold;

            // iterate backward until a non-tie is found
            while (true)
            {
                const double& i = thresholdItr->intensity;
                if (thresholdItr == mzIntensityPairs.begin())
                    break;
                else if (i != (--thresholdItr)->intensity)
                {
                    ++thresholdItr;
                    break;
                }
            }
            break;

        case ThresholdingBy_CountAfterTies:
            // no need to check bounds on thresholdItr because it gets checked above
            thresholdItr = mzIntensityPairs.begin() + ((size_t) threshold)-1;

            // iterate forward until a non-tie is found
            while (true)
            {
                const double& i = thresholdItr->intensity;
                if (++thresholdItr == mzIntensityPairs.end() ||
                    i != thresholdItr->intensity)
                    break;
            }
            break;

        case ThresholdingBy_AbsoluteIntensity:
            if (orientation == Orientation_MostIntense)
                thresholdItr = lower_bound(mzIntensityPairs.begin(),
                                           mzIntensityPairs.end(),
                                           MZIntensityPair(0, threshold),
                                           orientationMore_Predicate);
            else
                thresholdItr = lower_bound(mzIntensityPairs.begin(),
                                           mzIntensityPairs.end(),
                                           MZIntensityPair(0, threshold),
                                           orientationLess_Predicate);
            break;

        case ThresholdingBy_FractionOfBasePeakIntensity:
            if (orientation == Orientation_MostIntense)
                thresholdItr = lower_bound(mzIntensityPairs.begin(),
//.........这里部分代码省略.........

    void writeSpectrumText(SpectrumPtr s, ostream& os)
    {
        os << std::setprecision(7); // 123.4567
        bool ms1File = s->cvParam(MS_ms_level).valueAs<int>() == 1;
        
        // Write the scan numbers 
        os << "S\t";
        int scanNum = getScanNumber(s);
        os << scanNum <<  "\t" << scanNum;

        if (!ms1File)
        {
            // Write the precursor mz
            Precursor& precur = s->precursors[0];
            double mz = precur.isolationWindow.cvParam(MS_isolation_window_target_m_z).valueAs<double>();
            os << "\t" << mz;
        }
        os << "\n";
        
        // Write the scan time, if available
        if( !(s->scanList.empty()) && s->scanList.scans[0].cvParam(MS_scan_start_time).timeInSeconds() )
          os << "I\tRTime\t" << s->scanList.scans[0].cvParam(MS_scan_start_time).timeInSeconds()/60 << "\n";

        if (s->defaultArrayLength > 0)
        {
            // Write the base peak intensity and base peak m/z
            if (s->hasCVParam(MS_base_peak_intensity))
            {
                double bpi = s->cvParam(MS_base_peak_intensity).valueAs<double>();
                os << "I\tBPI\t" << bpi << "\n";
            }
            if (s->hasCVParam(MS_base_peak_m_z))
            {
                double bpm = s->cvParam(MS_base_peak_m_z).valueAs<double>();
                os << "I\tBPM\t" << bpm << "\n";
            }

            // Write the total ion current
            if (s->hasCVParam(MS_total_ion_current))
            {
                double tic = s->cvParam(MS_total_ion_current).valueAs<double>();
                os << "I\tTIC\t" << tic << "\n";
            }
        }

        //TODO 
        // Write ConvA/ConvB information
        // Write the ion injection time
        //os << "I\tIIT\t" << s->scanList.scans[0].cvParam(MS_ion_injection_time).timeInSeconds()/60 << "\n";

        if (!ms1File)
        {
            Precursor& precur = s->precursors[0];
            SelectedIon& si = precur.selectedIons[0];
            // Collect charge and mass info
            vector<int> charges;
            vector<double> masses;
            int numChargeStates = 0;
            // for each selected ion
            BOOST_FOREACH(const SelectedIon& curIon, precur.selectedIons){
              numChargeStates += getChargeStates(curIon, charges, masses);
            }

            // Write EZ lines if accurate masses are available
            CVParam massParam = si.cvParam(MS_accurate_mass_OBSOLETE);
            if( !massParam.empty() ){
              for(int i=0; i < numChargeStates; i++){
                os << "I\tEZ\t" << charges[i] << "\t" << masses[i] << "\t0\t0" << endl; // pad last two fields with 0
              }
            }

            // For each charge, write the charge and mass
            for(int i = 0; i < numChargeStates; i++)
            {
              os << "Z\t" << charges[i] << "\t" << masses[i] << "\n"; 
            }
        }

void test(bool indexed)
{
    if (os_) *os_ << "test(): indexed=\"" << boolalpha << indexed << "\"\n";

    MSData tiny;
    examples::initializeTiny(tiny);

    Serializer_mzML::Config config;
    config.indexed = indexed;
    Serializer_mzML serializer(config);  

    ostringstream oss;
    serializer.write(oss, tiny);

    if (os_) *os_ << "oss:\n" << oss.str() << endl;

    shared_ptr<istream> is(new istringstream(oss.str()));

    // dummy would normally be read in from file
  
    MSData dummy;

    ParamGroupPtr pg1(new ParamGroup);
    pg1->id = "CommonMS1SpectrumParams";
    pg1->cvParams.push_back(MS_positive_scan);
    pg1->cvParams.push_back(MS_full_scan);
    dummy.paramGroupPtrs.push_back(pg1);

    ParamGroupPtr pg2(new ParamGroup);
    pg2->id = "CommonMS2SpectrumParams";
    pg2->cvParams.push_back(MS_positive_scan);
    pg2->cvParams.push_back(MS_full_scan);
    dummy.paramGroupPtrs.push_back(pg2);

    // so we don't have any dangling references
    dummy.instrumentConfigurationPtrs.push_back(InstrumentConfigurationPtr(new InstrumentConfiguration("LCQDeca")));
    dummy.dataProcessingPtrs.push_back(DataProcessingPtr(new DataProcessing("XcaliburProcessing")));

    SpectrumListPtr sl = SpectrumList_mzML::create(is, dummy, indexed);

    // check easy functions

    unit_assert(sl.get());
    unit_assert(sl->size() == 4);
    unit_assert(sl->find ("S19") == 0);
    unit_assert(sl->findNative("19") == 0);
    unit_assert(sl->find("S20") == 1);
    unit_assert(sl->findNative("20") == 1);
    unit_assert(sl->find("S21") == 2);
    unit_assert(sl->findNative("21") == 2);
    unit_assert(sl->find("S22") == 3);
    unit_assert(sl->findNative("22") == 3);

    unit_assert(sl->findSpotID("A1").empty());
    IndexList spotIndexList = sl->findSpotID("A1,42x42,4242x4242");
    unit_assert(spotIndexList.size() == 1);
    unit_assert(spotIndexList[0] == 3);


    // check scan 19

    SpectrumPtr s = sl->spectrum(0); // read without binary data
    unit_assert(s.get());
    unit_assert(s->id == "S19");
    unit_assert(s->nativeID == "19");
    unit_assert(s->spotID.empty());
    unit_assert(s->cvParam(MS_ms_level).valueAs<int>() == 1);
    unit_assert(s->binaryDataArrayPtrs.empty());

    unit_assert(sl->spectrumIdentity(0).index == 0);
    unit_assert(sl->spectrumIdentity(0).id == "S19");
    unit_assert(sl->spectrumIdentity(0).nativeID == "19");
    unit_assert(sl->spectrumIdentity(0).spotID.empty());
 
    SpectrumPtr s_cache = sl->spectrum(0); // cache read
    unit_assert(s_cache.get() == s.get());

    s = sl->spectrum(0, true); // read with binary data
    unit_assert(s_cache.get() != s.get());

    vector<MZIntensityPair> pairs;
    s->getMZIntensityPairs(pairs);
    unit_assert(pairs.size() == 15);
    for (int i=0; i<15; i++)
        unit_assert(pairs[i].mz==i && pairs[i].intensity==15-i);

    unit_assert(s->spectrumDescription.scan.paramGroupPtrs.size() == 1);
    unit_assert(s->spectrumDescription.scan.paramGroupPtrs.back()->id == "CommonMS1SpectrumParams");
    unit_assert(s->spectrumDescription.scan.paramGroupPtrs.back()->cvParams.size() == 2);

    // check scan 20

    s = sl->spectrum(1, true);
    unit_assert(s.get());
    unit_assert(s->id == "S20");
    unit_assert(s->nativeID == "20");
    unit_assert(s->spotID.empty());
    unit_assert(s->cvParam(MS_ms_level).valueAs<int>() == 2);

    unit_assert(sl->spectrumIdentity(1).index == 1);
//.........这里部分代码省略.........