pmcsqs.cpp
来自「MS-Clustering is designed to rapidly clu」· C++ 代码 · 共 1,870 行 · 第 1/4 页
CPP
1,870 行
num_strong_pairs++;
inten_both+= peaks[i].intensity + peaks[didx].intensity;
}
}
sam.f_vals.push_back(fval(SQS_PROP_PEAKS_WITH_C1C2,num_pairs*one_over_np));
sam.f_vals.push_back(fval(SQS_PROP_STRONG_PEAKS_WITH_C1C2,num_strong_pairs*one_over_ns));
sam.f_vals.push_back(fval(SQS_PROP_INTEN_WITH_C1C2,inten_both*one_over_total_intensity));
}
// tag features
if (1)
{
vector<DPColumn> dp;
vector<float> tmp_vals;
int frag_charge;
for (frag_charge=1; frag_charge<=2; frag_charge++)
{
fill_SQS_DP(bs,dp,frag_charge);
float inten_in_tags[3]={0,0,0};
int num_peaks_in_tags[3]={0,0,0};
int num_strong_peaks_in_tags[3]={0,0,0};
vector<int> tag_lengths;
tag_lengths.resize(num_peaks,0);
int longest=0;
int i;
for (i=1; i<num_peaks; i++)
{
int max_tl=-1;
int aa;
for (aa=Ala; aa<=Val; aa++)
{
const int prev = dp[i].pointers[aa];
if (prev>=0 && tag_lengths[prev]>max_tl)
max_tl = tag_lengths[prev];
}
tag_lengths[i]=max_tl+1;
if (tag_lengths[i]>longest)
longest = tag_lengths[i];
int k;
for (k=0; k<3; k++)
{
if (tag_lengths[i]==k)
break;
inten_in_tags[k]+=peaks[i].intensity;
num_peaks_in_tags[k]++;
if (curr_spec_strong_inds[i])
num_strong_peaks_in_tags[k]++;
}
}
const int idx_off = (SQS_C2_IND_MAX_TAG_LENGTH_ABOVE_4 - SQS_IND_MAX_TAG_LENGTH_ABOVE_4)*(frag_charge-1);
if (longest>=4)
{
sam.f_vals.push_back(fval(SQS_IND_MAX_TAG_LENGTH_ABOVE_4+idx_off,1.0));
sam.f_vals.push_back(fval(SQS_MAX_TAG_LENGTH_ABOVE_4+idx_off,(float)longest-4.0));
}
else
{
sam.f_vals.push_back(fval(SQS_IND_MAX_TAG_LENGTH_BELOW_4+idx_off,1.0));
sam.f_vals.push_back(fval(SQS_MAX_TAG_LENGTH_BELOW_4+idx_off,(float)longest));
}
float inten_in_tags_both_sides=0;
for (i=0; i<num_peaks; i++)
if (dp[i].pointers[0])
inten_in_tags_both_sides+=peaks[i].intensity;
sam.f_vals.push_back(fval(SQS_PROP_INTEN_IN_TAGS+idx_off,inten_in_tags_both_sides*one_over_total_intensity));
tmp_vals.push_back(longest);
tmp_vals.push_back(inten_in_tags_both_sides*one_over_total_intensity);
// cout << "MAX TAG: " << longest << endl;
// cout << "Inten in TAG: " << inten_in_tags_both_sides/total_intensity << endl;
const float strong_threshes[]={0.3,0.2,0.1};
int k;
for (k=0; k<3; k++)
{
const int pos_off = 4*k;
const float prop_tags = (float)num_peaks_in_tags[k]*one_over_np;
const float prop_strong = (float)num_strong_peaks_in_tags[k]*one_over_ns;
const float prop_inten = (float)inten_in_tags[k]*one_over_total_intensity;
sam.f_vals.push_back(fval(SQS_PROP_TAGS1+idx_off+pos_off,prop_tags));
sam.f_vals.push_back(fval(SQS_PROP_STRONG_PEAKS_IN_TAG1+idx_off+pos_off,prop_strong));
sam.f_vals.push_back(fval(SQS_PROP_INTEN_TAG1+idx_off+pos_off,prop_inten));
if (prop_strong<strong_threshes[k])
sam.f_vals.push_back(fval(SQS_IND_PROP_STRONG_BELOW30_TAG1+idx_off+pos_off,1.0));
// save vals for diff features
tmp_vals.push_back(prop_tags);
tmp_vals.push_back(prop_strong);
tmp_vals.push_back(prop_inten);
}
}
int pos_off = tmp_vals.size()/2;
int i;
for (i=0; i<pos_off; i++)
sam.f_vals.push_back(fval(SQS_DIFF_MAX_TAG_LENGTH+i,tmp_vals[i]-tmp_vals[pos_off+i]));
/*
SQS_DIFF_MAX_TAG_LENGTH, SQS_DIFF_PROP_INTEN_IN_TAGS,
SQS_DIFF_PROP_TAGS1, SQS_DIFF_PROP_STRONG_PEAKS_IN_TAG1, SQS_DIFF_PROP_INTEN_TAG1,
SQS_DIFF_PROP_TAGS2, SQS_DIFF_PROP_STRONG_PEAKS_IN_TAG2, SQS_DIFF_PROP_INTEN_TAG2,
SQS_DIFF_PROP_TAGS3, SQS_DIFF_PROP_STRONG_PEAKS_IN_TAG3, SQS_DIFF_PROP_INTEN_TAG3,
*/
}
// density features
if (1)
{
mass_t t1=m_over_z*0.6666;
mass_t t2=m_over_z*1.3333;
mass_t h1=m_over_z;
float inten_t1=0,inten_t2=0, inten_h1=0;
int num_peaks_t1=0, num_peaks_t2=0, num_peaks_h1=0;
int i;
for (i=0; i<num_peaks && peaks[i].mass<t1; i++)
inten_t1+=peaks[i].intensity;
num_peaks_t1=i;
for ( ; i<num_peaks && peaks[i].mass<h1; i++)
inten_h1+=peaks[i].intensity;
num_peaks_h1=i;
inten_t2=inten_h1;
for ( ; i<num_peaks && peaks[i].mass<t2; i++)
inten_t2+=peaks[i].intensity;
num_peaks_t2=i-num_peaks_t1;
int num_peaks_t3= num_peaks - num_peaks_t2 - num_peaks_t1;
float inten_t3 = total_intensity - inten_t2 - inten_t1;
sam.f_vals.push_back(fval(SQS_PEAK_DENSE_T1,(float)(num_peaks_t1)*one_over_np));
sam.f_vals.push_back(fval(SQS_PEAK_DENSE_T2,(float)(num_peaks_t2)*one_over_np));
sam.f_vals.push_back(fval(SQS_PEAK_DENSE_T3,(float)(num_peaks_t3)*one_over_np));
// cout << "T1 PEAKS: " << (float)(num_peaks_t1)/(float)num_peaks << endl;
// cout << "T2 PEAKS: " << (float)(num_peaks_t2)/(float)num_peaks << endl;
// cout << "T3 PEAKS: " << (float)(num_peaks_t3)/(float)num_peaks << endl;
sam.f_vals.push_back(fval(SQS_INTEN_DENSE_T1,inten_t1*one_over_total_intensity));
sam.f_vals.push_back(fval(SQS_INTEN_DENSE_T2,inten_t2*one_over_total_intensity));
sam.f_vals.push_back(fval(SQS_INTEN_DENSE_T3,inten_t3*one_over_total_intensity));
// cout << "T1 INTEN: " << inten_t1/total_intensity << endl;
// cout << "T2 INTEN: " << inten_t2/total_intensity << endl;
// cout << "T3 INTEN: " << inten_t3/total_intensity << endl;
int num_peaks_h2 = num_peaks - num_peaks_h1;
float inten_h2 = total_intensity - inten_h1;
sam.f_vals.push_back(fval(SQS_PEAK_DENSE_H1,(num_peaks_h1)*one_over_np));
sam.f_vals.push_back(fval(SQS_PEAK_DENSE_H2,(num_peaks_h2)*one_over_np));
sam.f_vals.push_back(fval(SQS_INTEN_DENSE_H1,inten_h1*one_over_total_intensity));
sam.f_vals.push_back(fval(SQS_INTEN_DENSE_H2,inten_h2*one_over_total_intensity));
const float inten33=0.333* total_intensity;
const float inten50=0.5 * total_intensity;
const float inten75=0.75 * total_intensity;
const float inten90=0.90 * total_intensity;
float cum_inten=0;
const mass_t one_over_mz = 1.0 / (m_over_z + 0.1);
for (i=0; i<num_peaks && cum_inten<inten33; i++)
cum_inten+=peaks[i].intensity;
if (i==num_peaks)
i--;
sam.f_vals.push_back(fval(SQS_PROP_MZ_RANGE_WITH_33_INTEN,(peaks[i].mass*one_over_mz)-0.333));
for ( ; i<num_peaks && cum_inten<inten50; i++)
cum_inten+=peaks[i].intensity;
if (i==num_peaks)
i--;
sam.f_vals.push_back(fval(SQS_PROP_MZ_RANGE_WITH_50_INTEN,(peaks[i].mass*one_over_mz)-0.5));
for ( ; i<num_peaks && cum_inten<inten75; i++)
cum_inten+=peaks[i].intensity;
if (i==num_peaks)
i--;
sam.f_vals.push_back(fval(SQS_PROP_MZ_RANGE_WITH_75_INTEN,(peaks[i].mass*one_over_mz)-0.75));
for ( ; i<num_peaks && cum_inten<inten90; i++)
cum_inten+=peaks[i].intensity;
if (i==num_peaks)
i--;
sam.f_vals.push_back(fval(SQS_PROP_MZ_RANGE_WITH_90_INTEN,(peaks[i].mass*one_over_mz)-0.9));
}
// pmc features
if (1)
{
static vector< vector<float> > sqs_features;
get_sqs_features_from_pmc_tables(bs,sqs_features);
vector<float> max_vals;
max_vals.resize(4,0.000001);
vector< vector<float> > tmp_vals;
tmp_vals.resize(4);
// first give absolute counts
int charge;
for (charge=1; charge<=3; charge++)
{
vector<float>& vals = sqs_features[charge];
const int idx_off = 4 *( charge-1);
int i;
for (i=0; i<4; i++)
{
sam.f_vals.push_back(fval(SQS_NUM_FRAG_PAIRS_1+idx_off+i,vals[i]));
tmp_vals[charge].push_back(vals[i]);
if (vals[i]>max_vals[i])
max_vals[i]=vals[i];
}
}
// add features for prop of max
for (charge=1; charge<=3; charge++)
{
vector<float>& vals = sqs_features[charge];
const int idx_off = 4 *( charge-1);
int i;
for (i=0; i<4; i++)
{
const float ratio = (max_vals[i]>0 ? vals[i]/max_vals[i] : 0);
sam.f_vals.push_back(fval(SQS_PROP_OF_MAX_FRAG_PAIRS_1+idx_off+i,ratio));
tmp_vals[charge].push_back(ratio);
}
}
// conver to proportions by dividing by the number of peaks/strong_peaks
// and subtract the background levels (first)
const float one_over_peaks = 1.0/((float)num_peaks+0.1);
const float one_over_strong = 1.0/((float)num_strong+0.1);
for (charge=1; charge<=3; charge++)
{
vector<float>& vals = sqs_features[charge];
const int idx_off = 4 *( charge-1);
// normalize to get proportions of total number of peaks/strong
vals[0]-= this->curr_spec_rank_background_stats[charge].num_frag_pairs;
vals[1]-= this->curr_spec_rank_background_stats[charge].num_strong_frag_pairs;
vals[2]-= this->curr_spec_rank_background_stats[charge].num_c2_frag_pairs;
vals[3]-= this->curr_spec_rank_background_stats[charge].num_strong_c2_frag_pairs;
vals[0]*=one_over_peaks;
vals[1]*=one_over_strong;
vals[2]*=one_over_peaks;
vals[3]*=one_over_strong;
int i;
for (i=0; i<4; i++)
{
sam.f_vals.push_back(fval(SQS_PROP_FRAG_PAIRS_1+idx_off+i,vals[i]));
tmp_vals[charge].push_back(vals[i]);
}
}
// diff is between values for c2 and c3
int i;
for (i=0; i<tmp_vals[2].size(); i++)
sam.f_vals.push_back(fval(SQS_DIFF_NUM_FRAG_PAIRS_23+i,tmp_vals[2][i]-tmp_vals[3][i]));
}
sort(sam.f_vals.begin(),sam.f_vals.end());
}
/*******************************************************************************
Calculates the average statistics observed for the background using different
mass offsets
********************************************************************************/
void calc_background_stats(const mass_t single_charge_pair_sum, // the sum of b+y or c+z
Config *config,
const QCPeak *peaks,
const int num_peaks,
const vector<bool>& strong_inds,
const vector<float>& iso_levels,
const vector<bool>& strict_iso_inds,
PMCRankStats& pmc_stats_total)
{
const mass_t tolerance = config->get_tolerance();
const mass_t offsets[]={-22.0,-10.0,-8.5,8.5,12.0,22.5};
const int num_offsets = sizeof(offsets)/sizeof(mass_t);
const float one_over_num_offsets = 1.0 / (float)num_offsets;
int i;
for (i=0; i<num_offsets; i++)
{
static PMCRankStats pmc_stats;
calc_pmc_rank_stats_for_mass(peaks,num_peaks, single_charge_pair_sum+offsets[i],
1.5*tolerance,iso_levels,strong_inds, strict_iso_inds,pmc_stats);
if (i==0)
{
pmc_stats_total = pmc_stats;
}
else
{
pmc_stats_total.num_c2_frag_pairs += pmc_stats.num_c2_frag_pairs;
pmc_stats_total.num_frag_pairs += pmc_stats.num_frag_pairs;
pmc_stats_total.num_strong_c2_frag_pairs += pmc_stats.num_strong_c2_frag_pairs;
pmc_stats_total.num_strong_frag_pairs += pmc_stats.num_strong_frag_pairs;
pmc_stats_total.inten_frag_pairs += pmc_stats.inten_frag_pairs;
pmc_stats_total.inten_strong_pairs += pmc_stats.inten_strong_pairs;
pmc_stats_total.inten_c2_pairs += pmc_stats.inten_c2_pairs;
pmc_stats_total.inten_c2_strong_pairs += pmc_stats.inten_c2_strong_pairs;
}
}
pmc_stats_total.num_c2_frag_pairs *= one_over_num_offsets;
pmc_stats_total.num_frag_pairs *= one_over_num_offsets;
pmc_stats_total.num_strong_c2_frag_pairs *= one_over_num_offsets;
pmc_stats_total.num_strong_frag_pairs *= one_over_num_offsets;
pmc_stats_total.inten_frag_pairs *= one_over_num_offsets;
pmc_stats_total.inten_strong_pairs *= one_over_num_offsets;
pmc_stats_total.inten_c2_pairs *= one_over_num_offsets;
pmc_stats_total.inten_c2_strong_pairs *= one_over_num_offsets;
}
/*********************************************************************************
Fills the PMC sample features
Assumes the frag_pair_sum_offset is set.
Sets the PMC features and also creates vetors of features that are to be added
to the SQS sampels.
**********************************************************************************/
void PMCSQS_Scorer::calculate_curr_spec_pmc_values(const BasicSpectrum& bs,
mass_t increment)
{
const mass_t m_over_z = bs.ssf->m_over_z;
if (frag_pair_sum_offset<-999)
{
cout << "Error: must first set the expected frag pair offset!" << endl;
exit(1);
}
int charge;
for (charge=1; charge<=max_model_charge; charge++)
{
const mass_t org_pm_with_19 = bs.ssf->m_over_z * charge - 1.0025*(charge - 1);
const mass_t frag_pair_sum = org_pm_with_19 + frag_pair_sum_offset;
const int size_idx = this->get_rank_model_size_idx(charge,org_pm_with_19);
float bias = 0;
if (pmc_charge_mz_biases[charge].size()>size_idx)
bias = this->pmc_charge_mz_biases[charge][size_idx];
int bin_range =3*charge;
if (charge>=2)
bin_range = int(2*charge);
if (charge>=3)
bin_range = int(1.5*charge);
fill_rank_PMC_stats( charge,
frag_pair_sum + bias,
-bin_range-1,
bin_range-1,
increment,
config,
bs,
curr_spec_strong_inds,
curr_spec_iso_levels,
curr_spec_strict_iso_ind,
curr_spec_rank_pmc_tables[charge]);
calc_background_stats(frag_pair_sum,config,bs.peaks,bs.num_peaks,
curr_spec_strong_inds, curr_spec_iso_levels,
this->curr_spec_strict_iso_ind, curr_spec_rank_background_stats[charge]);
// find maximal values
int i;
PMCRankStats& maximal = curr_spec_rank_maximal_values[charge];
maximal = curr_spec_rank_pmc_tables[charge][0];
for (i=1; i<curr_spec_rank_pmc_tables[charge].size(); i++)
{
const PMCRankStats& curr_stats = curr_spec_rank_pmc_tables[charge][i];
// mathced pairs look for maximum number
if (curr_stats.num_c2_frag_pairs>maximal.num_c2_frag_pairs)
maximal.num_c2_frag_pairs = curr_stats.num_c2_frag_pairs;
if (curr_stats.num_frag_pairs>maximal.num_frag_pairs)
maximal.num_frag_pairs = curr_stats.num_frag_pairs;
if (curr_stats.num_strong_frag_pairs>maximal.num_strong_frag_pairs)
maximal.num_strong_frag_pairs = curr_stats.num_strong_frag_pairs;
if (curr_stats.num_strong_c2_frag_pairs > maximal.num_strong_c2_frag_pairs)
maximal.num_strong_c2_frag_pairs = curr_stats.num_strong_c2_frag_pairs;
// intensity look for maximal numbers
if (curr_stats.inten_frag_pairs > maximal.inten_frag_pairs)
maximal.inten_frag_pairs = curr_stats.inten_frag_pairs;
if (curr_stats.inten_strong_pairs> maximal.inten_strong_pairs)
maximal.inten_strong_pairs = curr_stats.inten_strong_pairs;
if (curr_stats.inten_c2_pairs> maximal.inten_c2_pairs)
maximal.inten_c2_pairs = curr_stats.inten_c2_pairs;
if (curr_stats.inten_c2_strong_pairs > maximal.inten_c2_strong_pairs)
maximal.inten_c2_strong_pairs = curr_stats.inten_c2_strong_pairs;
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?