C++ PeakMap::getPrimaryMSRunPath方法代码示例

  ExitCodes main_(int, const char **) override
  {
    //-------------------------------------------------------------
    // parameter handling
    //-------------------------------------------------------------

    //input/output files
    String in(getStringOption_("in"));
    String out(getStringOption_("out"));

    //-------------------------------------------------------------
    // loading input
    //-------------------------------------------------------------

    PeakMap exp;
    MzMLFile f;
    f.setLogType(log_type_);

    PeakFileOptions options;
    options.clearMSLevels();
    options.addMSLevel(2);
    f.getOptions() = options;
    f.load(in, exp);

    writeDebug_("Data set contains " + String(exp.size()) + " spectra", 1);

    //-------------------------------------------------------------
    // calculations
    //-------------------------------------------------------------

    vector<PeptideIdentification> pep_ids;
    CompNovoIdentificationCID comp_novo_id;

    // set the options
    Param algorithm_param = getParam_().copy("algorithm:", true);
    comp_novo_id.setParameters(algorithm_param);
    comp_novo_id.getIdentifications(pep_ids, exp);
    algorithm_param = comp_novo_id.getParameters();

    //-------------------------------------------------------------
    // writing output
    //-------------------------------------------------------------

    DateTime now = DateTime::now();
    String date_string = now.get();
    String identifier("CompNovoCID_" + date_string);

    for (vector<PeptideIdentification>::iterator it = pep_ids.begin(); it != pep_ids.end(); ++it)
    {
      it->assignRanks();
      it->setIdentifier(identifier);
    }

    vector<ProteinIdentification> prot_ids;
    ProteinIdentification prot_id;
    prot_id.setIdentifier(identifier);
    prot_id.setDateTime(now);
    StringList ms_runs;
    exp.getPrimaryMSRunPath(ms_runs);
    prot_id.setPrimaryMSRunPath(ms_runs);

    ProteinIdentification::SearchParameters search_parameters;
    search_parameters.charges = "+2-+3";
    if (algorithm_param.getValue("tryptic_only").toBool())
    {
      search_parameters.digestion_enzyme = *(ProteaseDB::getInstance()->getEnzyme("Trypsin"));
    }
    else
    {
      search_parameters.digestion_enzyme = *(ProteaseDB::getInstance()->getEnzyme("no cleavage"));
    }
    search_parameters.mass_type = ProteinIdentification::MONOISOTOPIC;
    search_parameters.fixed_modifications = algorithm_param.getValue("fixed_modifications");
    search_parameters.variable_modifications = algorithm_param.getValue("variable_modifications");

    search_parameters.missed_cleavages = (UInt)algorithm_param.getValue("missed_cleavages");
    search_parameters.fragment_mass_tolerance = (double)algorithm_param.getValue("fragment_mass_tolerance");
    search_parameters.precursor_mass_tolerance = (double)algorithm_param.getValue("precursor_mass_tolerance");
    search_parameters.fragment_mass_tolerance_ppm = false;
    search_parameters.precursor_mass_tolerance_ppm = false;
    prot_id.setSearchParameters(search_parameters);
    prot_id.setSearchEngineVersion("0.9beta");
    prot_id.setSearchEngine("CompNovo");
    prot_ids.push_back(prot_id);

    IdXMLFile().store(out, prot_ids, pep_ids);

    return EXECUTION_OK;
  }

//.........这里部分代码省略.........
    parameters << "-cfg" << cfg_file;

    // path to input file must be the last parameter
    parameters << inputfile_name;

    //-------------------------------------------------------------
    // calculations
    //-------------------------------------------------------------
    QStringList qparam;
    writeDebug_("MyriMatch arguments:", 1);
    writeDebug_(String("\"") + ListUtils::concatenate(parameters, "\" \"") + "\"", 1);
    for (Size i = 0; i < parameters.size(); ++i)
    {
      qparam << parameters[i].toQString();
    }
    QProcess process;

    // Bad style, because it breaks relative paths?
    process.setWorkingDirectory(tmp_dir.toQString());

    process.start(myrimatch_executable.toQString(), qparam, QIODevice::ReadOnly);
    bool success = process.waitForFinished(-1);
    String myri_msg(QString(process.readAllStandardOutput()));
    String myri_err(QString(process.readAllStandardError()));
    writeDebug_(myri_msg, 1);
    writeDebug_(myri_err, 0);
    if (!success || process.exitStatus() != 0 || process.exitCode() != 0)
    {
      writeLog_("Error: MyriMatch problem! (Details can be seen in the logfile: \"" + logfile + "\")");
      writeLog_("Note: This message can also be triggered if you run out of space in your tmp directory");
      return EXTERNAL_PROGRAM_ERROR;
    }

    //-------------------------------------------------------------
    // reading MyriMatch output
    //-------------------------------------------------------------

    writeDebug_("Reading output of MyriMatch", 5);
    String exp_name = File::basename(inputfile_name);
    String pep_file = tmp_dir + File::removeExtension(exp_name) + ".pepXML";

    vector<ProteinIdentification> protein_identifications;
    vector<PeptideIdentification> peptide_identifications;

    PeakMap exp;
    if (File::exists(pep_file))
    {
      MzMLFile fh;
      fh.load(inputfile_name, exp);

      SpectrumMetaDataLookup lookup;
      lookup.readSpectra(exp.getSpectra());
      PepXMLFile().load(pep_file, protein_identifications,
                        peptide_identifications, exp_name, lookup);
    }
    else
    {
      writeLog_("Error: MyriMatch problem! No pepXML output file (expected as '" + pep_file + "') was generated by MyriMatch.");
      writeLog_("Note: This message can be triggered if no MS2 spectra were found or no identifications were made.");
      writeLog_("      Myrimatch expects MS2 spectra in mzML files to contain the MSn tag. MSSpectrum with MS level 2 is not sufficient. You can use FileConverter to create such an mzML file by converting from mzML --> mzXML --> mzML.");
      return EXTERNAL_PROGRAM_ERROR;
    }

    if (debug_level_ == 0)
    {
      QFile(pep_file.toQString()).remove();
      QFile(cfg_file.toQString()).remove();
    }
    else
    {
      writeDebug_(String("Not removing '") + pep_file + "' for debugging purposes. Please delete manually!", 1);
      writeDebug_(String("Not removing '") + cfg_file + "' for debugging purposes. Please delete manually!", 1);
    }
    //-------------------------------------------------------------
    // writing results
    //-------------------------------------------------------------
    ProteinIdentification::SearchParameters search_parameters;
    search_parameters.db = getStringOption_("database");
    ProteinIdentification::PeakMassType mass_type = getFlag_("precursor_mass_tolerance_avg") == true ? ProteinIdentification::AVERAGE : ProteinIdentification::MONOISOTOPIC;
    search_parameters.mass_type = mass_type;
    search_parameters.fixed_modifications = getStringList_("fixed_modifications");
    search_parameters.variable_modifications = getStringList_("variable_modifications");
    search_parameters.missed_cleavages = getIntOption_("MaxMissedCleavages");
    search_parameters.fragment_mass_tolerance = getDoubleOption_("fragment_mass_tolerance");
    search_parameters.precursor_mass_tolerance = getDoubleOption_("precursor_mass_tolerance");
    search_parameters.precursor_mass_tolerance_ppm = getStringOption_("precursor_mass_tolerance_unit") == "ppm" ? true : false;
    search_parameters.fragment_mass_tolerance_ppm = getStringOption_("fragment_mass_tolerance_unit") == "ppm" ? true : false;
    protein_identifications[0].setSearchParameters(search_parameters);
    protein_identifications[0].setSearchEngineVersion(myrimatch_version);
    protein_identifications[0].setSearchEngine("MyriMatch");

    if (!protein_identifications.empty())
    {
      StringList ms_runs;
      exp.getPrimaryMSRunPath(ms_runs);
      protein_identifications[0].setPrimaryMSRunPath(ms_runs);
    }
    IdXMLFile().store(outputfile_name, protein_identifications, peptide_identifications);
    return EXECUTION_OK;
  }

//.........这里部分代码省略.........
        for (FeatureMap::const_iterator fit = feature_map_light.begin(); fit != feature_map_light.end(); ++fit)
        {
          all_features.push_back(*fit);
          light.push_back(MatchedFeature(*fit, idx));
        }
        for (FeatureMap::const_iterator fit = feature_map_heavy.begin(); fit != feature_map_heavy.end(); ++fit)
        {
          all_features.push_back(*fit);
          heavy.push_back(MatchedFeature(*fit, idx));
        }

        if (!heavy.empty() && !light.empty())
        {
          writeDebug_("Finding best feature pair out of " + String(light.size()) + " light and " + String(heavy.size()) + " heavy matching features.", 1);
          // now find "good" matches, means the pair with the smallest m/z deviation
          Feature best_light, best_heavy;
          double best_deviation(numeric_limits<double>::max());
          Size best_idx(pairs.size());
          for (vector<MatchedFeature>::const_iterator fit1 = light.begin(); fit1 != light.end(); ++fit1)
          {
            for (vector<MatchedFeature>::const_iterator fit2 = heavy.begin(); fit2 != heavy.end(); ++fit2)
            {
              if (fit1->idx != fit2->idx || fit1->f.getCharge() != fit2->f.getCharge() ||
                  fabs(fit1->f.getMZ() - pairs[fit1->idx].mz_light) > precursor_mass_tolerance ||
                  fabs(fit2->f.getMZ() - pairs[fit2->idx].mz_heavy) > precursor_mass_tolerance)
              {
                continue;
              }
              double deviation(0);
              deviation = fabs((fit1->f.getMZ() - pairs[fit1->idx].mz_light) - (fit2->f.getMZ() - pairs[fit2->idx].mz_heavy));
              if (deviation < best_deviation && deviation < precursor_mass_tolerance)
              {
                best_light = fit1->f;
                best_heavy = fit2->f;
                best_idx = fit1->idx;
              }
            }
          }

          if (best_idx == pairs.size())
          {
            continue;
          }

          writeDebug_("Ratio: " + String(best_heavy.getIntensity() / best_light.getIntensity()), 1);
          ConsensusFeature SILAC_feature;
          SILAC_feature.setMZ((best_light.getMZ() + best_heavy.getMZ()) / 2.0);
          SILAC_feature.setRT((best_light.getRT() + best_heavy.getRT()) / 2.0);
          SILAC_feature.insert(0, best_light);
          SILAC_feature.insert(1, best_heavy);
          results_map.push_back(SILAC_feature);
          quantlets.push_back(SILACQuantitation(best_light.getIntensity(), best_heavy.getIntensity(), best_idx));
        }
      }
    }

    // now calculate the final quantitation values from the quantlets
    Map<Size, vector<SILACQuantitation> > idx_to_quantlet;
    for (vector<SILACQuantitation>::const_iterator it = quantlets.begin(); it != quantlets.end(); ++it)
    {
      idx_to_quantlet[it->idx].push_back(*it);
    }

    for (Map<Size, vector<SILACQuantitation> >::ConstIterator it1 = idx_to_quantlet.begin(); it1 != idx_to_quantlet.end(); ++it1)
    {
      SILAC_pair silac_pair = pairs[it1->first];

      // simply add up all intensities and calculate the final ratio
      double light_sum(0), heavy_sum(0);
      vector<double> light_ints, heavy_ints, ratios;
      for (vector<SILACQuantitation>::const_iterator it2 = it1->second.begin(); it2 != it1->second.end(); ++it2)
      {
        light_sum += it2->light_intensity;
        light_ints.push_back(it2->light_intensity);
        heavy_sum += it2->heavy_intensity;
        heavy_ints.push_back(it2->heavy_intensity);
        ratios.push_back(it2->heavy_intensity / it2->light_intensity * (it2->heavy_intensity + it2->light_intensity));
      }

      double absdev_ratios = Math::absdev(ratios.begin(), ratios.begin() + (ratios.size()) / (heavy_sum + light_sum));
      cout << "Ratio: " << silac_pair.mz_light << " <-> " << silac_pair.mz_heavy << " @ " << silac_pair.rt << " s, ratio(h/l) " << heavy_sum / light_sum << " +/- " << absdev_ratios << " (#scans for quantation: " << String(it1->second.size()) << " )" << endl;
    }


    //-------------------------------------------------------------
    // writing output
    //-------------------------------------------------------------
    StringList ms_runs;
    exp.getPrimaryMSRunPath(ms_runs);
    if (feature_out != "")
    {
      all_features.setPrimaryMSRunPath(ms_runs);
      FeatureXMLFile().store(feature_out, all_features);
    }
    writeDebug_("Writing output", 1);
    results_map.setPrimaryMSRunPath(ms_runs);
    ConsensusXMLFile().store(out, results_map);

    return EXECUTION_OK;
  }

//.........这里部分代码省略.........
            for (QStringList::ConstIterator subdir_file_it=subdir_files.begin(); subdir_file_it!=subdir_files.end(); ++subdir_file_it)
            {
              QFile::copy(qdir_models_source.filePath(*subdir_file_it), qdir_temp.filePath(*subdir_file_it));
            }
            qdir_temp.cdUp();
            qdir_models_source.cdUp();
          }
          else
          {
            QFile::copy(qdir_models_source.filePath(*file_it), qdir_temp.filePath(*file_it));
          }
        }

        //generate PTM File and store in temp directory
        PepNovoInfile p_novo_infile;
        String ptm_command;
        if (!getStringList_("fixed_modifications").empty() || !getStringList_("variable_modifications").empty())
        {
          p_novo_infile.setModifications(getStringList_("fixed_modifications"), getStringList_("variable_modifications"));
          p_novo_infile.store(qdir_temp.filePath("PepNovo_PTMs.txt"));
          pepnovo_files.append("PepNovo_PTMs.txt");
          p_novo_infile.getModifications(mods_and_keys);

          for (std::map<String, String>::const_iterator key_it=mods_and_keys.begin(); key_it!=mods_and_keys.end();++key_it)
          {
            if (ptm_command!="")
            {
              ptm_command+=":";
            }
            ptm_command+= key_it->first;
            //key_to_id[key_it->second]=key_it->first;
          }
        }

        //-------------------------------------------------------------
        // (3) running program according to parameters
        //-------------------------------------------------------------
        QStringList arguments;

        arguments << "-file" << mgf_file.toQString();
        arguments << "-model" << model_name.toQString();
        if (pm_tolerance != -1 ) arguments << "-pm_tolerance"<<String(pm_tolerance).toQString();
        if (fragment_tolerance != -1 ) arguments << "-fragment_tolerance" <<String(fragment_tolerance).toQString();
        if (!ptm_command.empty()) arguments <<"-PTMs" <<ptm_command.toQString();
        if (getFlag_("correct_pm")) arguments << "-correct_pm";
        if (getFlag_("use_spectrum_charge")) arguments << "-use_spectrum_charge";
        if (getFlag_("use_spectrum_mz")) arguments << "-use_spectrum_mz";
        if (getFlag_("no_quality_filter")) arguments << "-no_quality_filter";
        arguments << "-digest" << digest.toQString();
        arguments << "-num_solutions" << String(num_solutions).toQString();
        if (tag_length!=-1) arguments<<"-tag_length" << String(tag_length).toQString();
        arguments<<"-model_dir" << tmp_models_dir.toQString();
        //arguments<<">" << temp_pepnovo_outfile.toQString();

        writeDebug_("Use this line to call PepNovo: ", 1);
        writeDebug_(pepnovo_executable + " " + String(arguments.join(" ")), 1);
        QProcess process;
        process.setStandardOutputFile(temp_pepnovo_outfile.toQString());
        process.setStandardErrorFile(temp_pepnovo_outfile.toQString());
        process.start(pepnovo_executable.toQString(), arguments); // does automatic escaping etc...
        if (process.waitForFinished(-1))
        {
          //if PepNovo finished successfully use PepNovoOutfile to parse the results and generate idXML
          std::vector< PeptideIdentification > peptide_identifications;
          ProteinIdentification protein_identification;
          StringList ms_runs;
          exp.getPrimaryMSRunPath(ms_runs);
          protein_identification.setPrimaryMSRunPath(ms_runs);

          PepNovoOutfile p_novo_outfile;

          //resolve PTMs (match them back to the OpenMs Identifier String)
          std::vector<ProteinIdentification>prot_ids;
          p_novo_outfile.load(temp_pepnovo_outfile, peptide_identifications, protein_identification, -1e5, index_to_precursor, mods_and_keys);
          prot_ids.push_back(protein_identification);
          IdXMLFile().store(outputfile_name, prot_ids, peptide_identifications);
        }

        if (process.exitStatus() != 0)  error = true;
       
      }
      catch(Exception::BaseException &exc)
      {
        writeLog_(exc.what());
        LOG_ERROR << "Error occurred: " << exc.what() << std::endl;
        error = true;
      }
      
      if (!error)
      {
        File::removeDirRecursively(temp_data_directory);
        return EXECUTION_OK;
      }
      else
      {
        writeLog_("PepNovo problem. Aborting! (Details can be seen in outfiles: '" + temp_data_directory + "')");
        return EXTERNAL_PROGRAM_ERROR;
      }

    }

  ExitCodes main_(int, const char**) override
  {
    //input file names
    String in = getStringOption_("in");
    String out = getStringOption_("out");
    String out_mzq = getStringOption_("out_mzq");

    //prevent loading of fragment spectra
    PeakFileOptions options;
    options.setMSLevels(vector<Int>(1, 1));

    //reading input data
    MzMLFile f;
    f.getOptions() = options;
    f.setLogType(log_type_);

    PeakMap exp;
    f.load(in, exp);
    exp.updateRanges();

    if (exp.getSpectra().empty())
    {
      throw OpenMS::Exception::FileEmpty(__FILE__, __LINE__, __FUNCTION__, "Error: No MS1 spectra in input file.");
    }

    // determine type of spectral data (profile or centroided)
    SpectrumSettings::SpectrumType  spectrum_type = exp[0].getType();

    if (spectrum_type == SpectrumSettings::PROFILE)
    {
      if (!getFlag_("force"))
      {
        throw OpenMS::Exception::IllegalArgument(__FILE__, __LINE__, __FUNCTION__, "Error: Profile data provided but centroided spectra expected. To enforce processing of the data set the -force flag.");
      }
    }

    //load seeds
    FeatureMap seeds;
    if (getStringOption_("seeds") != "")
    {
      FeatureXMLFile().load(getStringOption_("seeds"), seeds);
    }

    //setup of FeatureFinder
    FeatureFinder ff;
    ff.setLogType(log_type_);

    // A map for the resulting features
    FeatureMap features;
    StringList ms_runs;
    exp.getPrimaryMSRunPath(ms_runs);
    features.setPrimaryMSRunPath(ms_runs);

    // get parameters specific for the feature finder
    Param feafi_param = getParam_().copy("algorithm:", true);
    writeDebug_("Parameters passed to FeatureFinder", feafi_param, 3);

    // Apply the feature finder
    ff.run(FeatureFinderAlgorithmPicked::getProductName(), exp, features, feafi_param, seeds);
    features.applyMemberFunction(&UniqueIdInterface::setUniqueId);

    // DEBUG
    if (debug_level_ > 10)
    {
      FeatureMap::Iterator it;
      for (it = features.begin(); it != features.end(); ++it)
      {
        if (!it->isMetaEmpty())
        {
          vector<String> keys;
          it->getKeys(keys);
          LOG_INFO << "Feature " << it->getUniqueId() << endl;
          for (Size i = 0; i < keys.size(); i++)
          {
            LOG_INFO << "  " << keys[i] << " = " << it->getMetaValue(keys[i]) << endl;
          }
        }
      }
    }

    //-------------------------------------------------------------
    // writing files
    //-------------------------------------------------------------

    //annotate output with data processing info
    addDataProcessing_(features, getProcessingInfo_(DataProcessing::QUANTITATION));

    // write features to user specified output file
    FeatureXMLFile map_file;

    // Remove detailed convex hull information and subordinate features
    // (unless requested otherwise) to reduce file size of feature files
    // unless debugging is turned on.
    if (debug_level_ < 5)
    {
      FeatureMap::Iterator it;
      for (it = features.begin(); it != features.end(); ++it)
      {
        it->getConvexHull().expandToBoundingBox();
        for (Size i = 0; i < it->getConvexHulls().size(); ++i)
//.........这里部分代码省略.........