本文整理汇总了C++中PeakMap::areaBeginConst方法的典型用法代码示例。如果您正苦于以下问题:C++ PeakMap::areaBeginConst方法的具体用法?C++ PeakMap::areaBeginConst怎么用?C++ PeakMap::areaBeginConst使用的例子?那么, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类PeakMap
的用法示例。
在下文中一共展示了PeakMap::areaBeginConst方法的1个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: filterByFoldChange
void filterByFoldChange(const PeakMap& exp1, const PeakMap& exp2,
const vector<double>& pc_ms2_rts, const vector<double>& pc_mzs,
const double rttol, const double mztol, double fold_change,
vector<double>& control_XIC_larger,
vector<double>& treatment_XIC_larger,
vector<double>& indifferent_XICs)
{
assert(pc_mzs.size() == pc_ms2_rts.size());
// search for each EIC and add up
for (Size i = 0; i < pc_mzs.size(); ++i)
{
//cerr << "start" << endl;
double pc_ms2_rt = pc_ms2_rts[i];
double pc_mz = pc_mzs[i];
//std::cerr << "Rt" << cm[i].getRT() << " mz: " << cm[i].getMZ() << " R " << cm[i].getMetaValue("rank") << "\n";
double mz_da = mztol * pc_mzs[i] / 1e6; // mz tolerance in Dalton
double rt_start = pc_ms2_rts[i] - rttol / 2.0;
// get area iterator (is MS1 only!) for rt and mz window
PeakMap::ConstAreaIterator it1 = exp1.areaBeginConst(pc_ms2_rt - rttol / 2, pc_ms2_rt + rttol / 2, pc_mz - mz_da, pc_mz + mz_da);
PeakMap::ConstAreaIterator it2 = exp2.areaBeginConst(pc_ms2_rt - rttol / 2, pc_ms2_rt + rttol / 2, pc_mz - mz_da, pc_mz + mz_da);
// determine maximum number of MS1 scans in retention time window
set<double> rts1;
set<double> rts2;
for (; it1 != exp1.areaEndConst(); ++it1)
{
rts1.insert(it1.getRT());
}
for (; it2 != exp2.areaEndConst(); ++it2)
{
rts2.insert(it2.getRT());
}
Size length = std::max(rts1.size(), rts2.size()) / 2.0;
//cout << length << endl;
if (length == 0)
{
cerr << "WARNING: no MS1 scans in retention time window found in both maps (mz: " << pc_mzs[i] << " / rt: " << pc_ms2_rts[i] << ")" << endl;
continue;
}
vector<double> XIC1(length, 0.0);
vector<double> XIC2(length, 0.0);
it1 = exp1.areaBeginConst(pc_ms2_rt - rttol / 2, pc_ms2_rt + rttol / 2, pc_mz - mz_da, pc_mz + mz_da);
it2 = exp2.areaBeginConst(pc_ms2_rt - rttol / 2, pc_ms2_rt + rttol / 2, pc_mz - mz_da, pc_mz + mz_da);
for (; it1 != exp1.areaEndConst(); ++it1)
{
double relative_rt = (it1.getRT() - rt_start) / rttol;
Size bin = relative_rt * (length - 1);
XIC1[bin] += it1->getIntensity();
if (bin >= length)
{
bin = length - 1;
}
}
for (; it2 != exp2.areaEndConst(); ++it2)
{
double relative_rt = (it2.getRT() - rt_start) / rttol;
Size bin = relative_rt * (length - 1);
if (bin >= length)
{
bin = length - 1;
}
XIC2[bin] += it2->getIntensity();
}
double total_itensity1 = std::accumulate(XIC1.begin(), XIC1.end(), 0.0);
double total_itensity2 = std::accumulate(XIC2.begin(), XIC2.end(), 0.0);
double ratio = total_itensity2 / (total_itensity1 + 1);
//cout << pc_ms2_rt << "/" << pc_mz << " has ratio: " << ratio << " determined on " << length << " bins" << endl;
if (ratio < 1.0 / fold_change)
{
control_XIC_larger.push_back(pc_ms2_rt);
}
else if (ratio > fold_change)
{
treatment_XIC_larger.push_back(pc_ms2_rt);
}
else
{
indifferent_XICs.push_back(pc_ms2_rt);
continue;
}
/*
for (Size k = 0; k != length; ++k)
{
cout << k << ": " << rt_start + rttol / length * k << ": " << XIC1[k] << " " << XIC2[k] << endl;
//.........这里部分代码省略.........