pmcsqs.cpp

来自「MS-Clustering is designed to rapidly clu」· C++ 代码 · 共 1,870 行 · 第 1/4 页

	
			if (j==pms_with_19.size())
			{
				pms_with_19.push_back(pm_with_19);
				charges.push_back(max_charge);

			//	cout << "1Add : " << pm_with_19 << " " << max_charge << endl;
			}
		}

		// add tol to -1 of 2nd
		if (charges[0]>=2)
		{
			const mass_t pm_with_19 = pms_with_19[1] - MASS_PROTON;
			int j;
			for (j=0; j<pms_with_19.size(); j++)
				if (charges[j]==max_charge &&
					fabs(pm_with_19 - pms_with_19[j])<half_tol)
						break;
			
			if (j==pms_with_19.size())
			{
				pms_with_19.push_back(pm_with_19);
				charges.push_back(max_charge);

			//	cout << "2Add : " << pm_with_19 << " " << max_charge << endl;
			}
		}
		

		if (charges[0]>=3)
		{
			const mass_t pm_with_19 = pms_with_19[1] + MASS_PROTON;
			int j;
			for (j=0; j<pms_with_19.size(); j++)
				if (charges[j]==max_charge &&
					fabs(pm_with_19 - pms_with_19[j])<half_tol)
						break;

			if (j==pms_with_19.size())
			{
				pms_with_19.push_back(pm_with_19);
				charges.push_back(max_charge);

			//	cout << "3Add : " << pm_with_19 << " " << max_charge << endl;
			}
		}
	}

	// add other charges if their comp probability and sqs probs are high enough
	if (specific_charge==0)
	{
		// find best charge
		int c; 
		float max_prob=-1.0;

		for (c=1; c<res.size(); c++)
		{
			if (c==max_charge)
				continue;

			if (res[c].min_comp_prob > min_comp_prob_for_adding_second &&
				res[c].sqs_prob > min_sqs_prob_for_adding_second)
			{
				pms_with_19.push_back(res[c].mz1 * c - (c-1)*MASS_PROTON);
				charges.push_back(c);

			//	cout << "4Add : " << res[c].mz1 * c - (c-1)*MASS_PROTON << " " << c << endl;
			}
		}
	}

	if (all_res)
		*all_res = res;
}



void PMCSQS_Scorer::benchmark_pm_selection(Config *config, FileManager& fm, mass_t pm_val_tol)
{
	const vector< vector< mass_t > >& threshes = config->get_size_thresholds();
	int c=1;

	for (c=1; c<threshes.size(); c++)
	{
		int s;
		for (s=0; s<threshes[c].size(); s++)
		{
			mass_t min_mz = (s>0 ? threshes[c][s-1]/c : 0);
			mass_t max_mz = threshes[c][s]/c;

			FileSet fs;
			fs.select_files_in_mz_range(fm,min_mz,max_mz,c);

			if (fs.get_total_spectra()<200)
				continue;

			cout << "CHARGE " << c <<" size " << s << "  (" << fs.get_total_spectra() << " spectra)" << endl;

			const vector<SingleSpectrumFile *>& all_ssfs = fs.get_ssf_pointers();
			BasicSpecReader bsr;
			vector<QCPeak> peaks;
			peaks.resize(10000);

			vector<int> correct_counts;
			correct_counts.resize(8,0);

			int num_correct=0;
			int num_wrong_charge=0;
			int num_diff_charge=0;

			int i;
			for (i=0; i<all_ssfs.size(); i++)
			{
				SingleSpectrumFile *ssf=all_ssfs[i];
				BasicSpectrum bs;

				bs.num_peaks = bsr.read_basic_spec(config,fm,ssf,&peaks[0]);
				bs.peaks = &peaks[0];
				bs.ssf = ssf;

				const mass_t true_mass = ssf->peptide.get_mass_with_19();

				vector<mass_t> pms_with_19;
				vector<int>    charges;
				select_pms_and_charges(config,bs,pms_with_19,charges,true);

				if (charges[0] != c)
					num_wrong_charge++;

				bool got_diff_charge=false;
				int j;
				for (j=1; j<charges.size(); j++)
					if (charges[j] != charges[0])
						got_diff_charge=true;

				if (got_diff_charge)
					num_diff_charge++;

				for (j=0; j<pms_with_19.size(); j++)
					if (fabs(pms_with_19[j]-true_mass)<pm_val_tol)
						break;
				
				if (j==pms_with_19.size())
					continue;
				

				num_correct++;

				if (got_diff_charge && j == pms_with_19.size()-1)
				{
					correct_counts[7]++;
				}
				else
					correct_counts[j]++;
			}

			double num_total = (double)all_ssfs.size();
			cout << "Had correct       " << fixed << setprecision(4) << num_correct/num_total << endl;
			cout << "First correct     " << correct_counts[0]/num_total << endl;
			cout << "Second correct    " << correct_counts[1]/num_total << endl;
			cout << "Off-1  correct    " << correct_counts[2]/num_total << endl;
			cout << "Off+1  correct    " << correct_counts[3]/num_total << endl;
			cout << "Off-1  2nd        " << correct_counts[4]/num_total << endl;
			cout << "Off+1  2nd        " << correct_counts[5]/num_total << endl;
			cout << "Diff Ch correct   " << correct_counts[7]/num_total << endl;
			cout << "With wrong charge " << num_wrong_charge/num_total << endl;
			cout << "With diff charge  " << num_diff_charge/num_total << endl << endl;
		}
	}
}


const int DPColumnBytes = sizeof(int)*(Val+1);

// for each peak,aa holds the idx of the previous aa if they have a mass
// diff of that aa (within tolerance)
// entry 0 in each column holds an indicator if peak is in aa diff
// entry 1 in each column holds an indicator if peak has a 
void PMCSQS_Scorer::fill_SQS_DP(const BasicSpectrum& bs, vector<DPColumn>& dp, int frag_charge ) const
{
	const QCPeak *peaks = bs.peaks;
	const int num_peaks = bs.num_peaks;
	const vector<mass_t>& aa2mass = config->get_aa2mass();
	const mass_t tag_tolerance = (config->get_tolerance()<0.15 ? 
									config->get_tolerance() : config->get_tolerance()*0.33);

	
	const mass_t mult_val = (1.0 / frag_charge);
	dp.resize(num_peaks);
	int i;
	for (i=0; i<num_peaks; i++)
	{
		dp[i].pointers[0]=0;
		int j;
		for (j=1; j<=Val; j++)
			dp[i].pointers[j]=-1;
	}

	int aa;
	for (aa=Ala; aa<=Val; aa++)
	{
		if (aa==Ile || aa==Xle)
			continue;

		const mass_t aamass = mult_val * aa2mass[aa];
		const mass_t min_offset = aamass-tag_tolerance;
		const mass_t max_offset = aamass+tag_tolerance;

		int trail_idx=0;
		int lead_idx=1;

		while (lead_idx<num_peaks)
		{
			if (curr_spec_iso_levels[lead_idx]>0)
			{
				lead_idx++;
				continue;
			}

			while (peaks[lead_idx].mass-peaks[trail_idx].mass>max_offset)
				trail_idx++;

			if (curr_spec_iso_levels[trail_idx]==0 && 
				peaks[lead_idx].mass-peaks[trail_idx].mass>min_offset)
			{
				dp[lead_idx].pointers[aa]=trail_idx;
				dp[lead_idx].pointers[0]=1;
				dp[trail_idx].pointers[0]=1;

			//	cout << "Off: " << peaks[lead_idx].mass-peaks[trail_idx].mass - min_offset<< endl;
			}
			else 
				dp[lead_idx].pointers[aa]=-1;

			lead_idx++;
		}

	}
}


/*****************************************************************************
Does the raw calculations required for processing the spectrum.
******************************************************************************/
bool PMCSQS_Scorer::init_for_current_spec(Config *_config, 
										  const BasicSpectrum& bs)
{
	config = _config;

	bs.calc_peak_isotope_levels(config->get_tolerance(),this->curr_spec_iso_levels);
	if (! bs.select_strong_peak_idxs(this->curr_spec_iso_levels,this->curr_spec_strong_inds))
		return false;

	bs.mark_all_possible_isotope_peaks(config->get_tolerance(),curr_spec_strict_iso_ind);

	curr_spec_total_intensity=0;
	curr_spec_strong_intensity=0;
	curr_spec_num_strong=0;
	int i;
	for (i=0; i<bs.num_peaks; i++)
	{
		if (curr_spec_iso_levels[i]>0)
			continue;

		curr_spec_total_intensity+=bs.peaks[i].intensity;
		if (curr_spec_strong_inds[i])
		{
			curr_spec_strong_intensity+=bs.peaks[i].intensity;
			curr_spec_num_strong++;
		}
	}


	if (curr_spec_rank_pmc_tables.size()<=max_model_charge+1)
	{
		curr_spec_rank_pmc_tables.clear();
		curr_spec_rank_pmc_tables.resize(max_model_charge+1);

		curr_spec_rank_background_stats.clear();
		curr_spec_rank_background_stats.resize(max_model_charge+1);

		curr_spec_rank_maximal_values.clear();
		curr_spec_rank_maximal_values.resize(max_model_charge+1);
	}

	return true;
}


/****************************************************************************
*****************************************************************************/
void PMCSQS_Scorer::fill_fval_vector_with_SQS(const BasicSpectrum& bs,
								   ME_Regression_Sample& sam) const
{
	const mass_t tolerance = config->get_tolerance();
	const mass_t frag_tolerance = (tolerance<0.15 ? tolerance : tolerance * 0.5);
	const QCPeak *peaks  = bs.peaks;
	const int  num_peaks = bs.num_peaks;
	const mass_t m_over_z = bs.ssf->m_over_z;
	const int num_strong=curr_spec_num_strong;
	const float total_intensity=curr_spec_total_intensity;
	const float one_over_np = 10.0 / (num_peaks+1);
	const float one_over_ns = 5.0/ (num_strong + 1);
	const float one_over_total_intensity = 1.0 / (total_intensity + 1.0);
	const mass_t max_peak_mass = (num_peaks <5 ? POS_INF : peaks[num_peaks-1].mass);

	sam.f_vals.clear();
	sam.f_vals.push_back(fval(SQS_CONST,1.0));
	sam.f_vals.push_back(fval(SQS_PEAK_DENSITY, (float)num_peaks/max_peak_mass));

	float grass_level_inten=NEG_INF;
	
	// calculate grass level peaks
	if (1)
	{
		int i;
		vector<float> peak_intens;
		peak_intens.resize(num_peaks);
		for (i=0; i<num_peaks; i++)
			peak_intens[i]=peaks[i].intensity;

		sort(peak_intens.begin(),peak_intens.end());

		int idx_G = num_peaks/3;
		grass_level_inten = peak_intens[idx_G];

		float cum_intensity2G=0;
		float inten_2G=2.0 * grass_level_inten;
		int idx_2G = idx_G;
		while (idx_2G<num_peaks && peak_intens[idx_2G]<inten_2G)
		{
			cum_intensity2G+=peak_intens[idx_2G];
			idx_2G++;
		}

		float cum_intensity5G=0;
		float inten_5G = 5.0 * grass_level_inten;
		int idx_5G = idx_2G;
		while (idx_5G<num_peaks && peak_intens[idx_5G]<inten_5G)
		{
			cum_intensity5G+=peak_intens[idx_5G];
			idx_5G++;
		}

		float inten_10G = 10.0 * grass_level_inten;
		int idx_10G = idx_5G;
		while (idx_10G<num_peaks && peak_intens[idx_10G]<inten_10G)
			idx_10G++;

		sam.f_vals.push_back(fval(SQS_PROP_UPTO2G,(float)idx_2G*one_over_np));
		sam.f_vals.push_back(fval(SQS_PROP_UPTO5G,(float)(idx_5G-idx_2G)*one_over_np));
		sam.f_vals.push_back(fval(SQS_PROP_UPTO10G,(float)(idx_10G-idx_5G)*one_over_np)); 
		sam.f_vals.push_back(fval(SQS_PROP_MORE10G,(float)(num_peaks-idx_10G)*one_over_np));
		sam.f_vals.push_back(fval(SQS_PROP_INTEN_UPTO2G,cum_intensity2G*one_over_total_intensity));
		sam.f_vals.push_back(fval(SQS_PROP_INTEN_UPTO5G,cum_intensity5G*one_over_total_intensity));
		sam.f_vals.push_back(fval(SQS_PROP_INTEN_MORE5G, (total_intensity-cum_intensity2G-cum_intensity5G)*one_over_total_intensity));
	}

	// isotope features
	if (1)
	{
		int i;
		int num_with_iso=0;
		int strong_with_iso=0;
		for (i=1; i<num_peaks; i++)
			if (curr_spec_iso_levels[i]>0)
			{
				num_with_iso++;
				if (curr_spec_iso_levels[i-1]==0 && curr_spec_strong_inds[i-1])
					strong_with_iso++;
			}

		sam.f_vals.push_back(fval(SQS_PROP_ISO_PEAKS,num_with_iso*one_over_np)); 
		sam.f_vals.push_back(fval(SQS_PROP_STRONG_WITH_ISO_PEAKS,strong_with_iso*one_over_ns)); 

	//	cout << "PROP WITH ISO  : " << num_with_iso/(float)num_peaks << endl;
	//	cout << "STRONG WITH IOS: " << strong_with_iso/(float)num_strong << endl;
	}


	// neutral loss features
	if (1)
	{
		vector<float> tmp_vals;
		tmp_vals.clear();
		int frag_charge;
		for (frag_charge=1; frag_charge<=2; frag_charge++)
		{
			const mass_t offsets[3]={MASS_H2O/frag_charge, MASS_NH3/frag_charge, MASS_CO/frag_charge};
			const int num_offsets = 3;
			
			const int fc_off = (frag_charge-1)*num_offsets * 2;
			int i;
			for (i=0; i<num_offsets; i++)
			{
				const mass_t min_offset = offsets[i]-frag_tolerance;
				const mass_t max_offset = offsets[i]+frag_tolerance;

				int num_pairs=0;
				int num_strong_pairs=0;

				int trail_idx=0;
				int lead_idx=1;

				while (lead_idx<num_peaks)
				{
					while (peaks[lead_idx].mass-peaks[trail_idx].mass>max_offset)
						trail_idx++;

					if (peaks[lead_idx].mass-peaks[trail_idx].mass>min_offset)
					{
						num_pairs++;
						if (curr_spec_strong_inds[lead_idx])
							num_strong_pairs++;
					}
					lead_idx++;
				}

				const float prop_peaks = num_pairs*one_over_np;
				const float prop_strong = num_strong_pairs*one_over_ns;
				sam.f_vals.push_back(fval(SQS_PROP_ALL_WITH_H2O_LOSS+fc_off+i,prop_peaks));
				sam.f_vals.push_back(fval(SQS_PROP_STRONG_WITH_H2O_LOSS+fc_off+i,prop_strong));

				tmp_vals.push_back(prop_peaks);
				tmp_vals.push_back(prop_strong);

			//	cout << "OFF REG " << i << " " << num_pairs/(float)num_peaks << endl;
			//	cout << "OFF STR " << i << " " << num_strong_pairs/(float)num_strong << endl;
			}
		}

		int i;
		const int half_size = tmp_vals.size()/2; 
		for (i=0; i<half_size; i++)
			sam.f_vals.push_back(fval(SQS_DIFF_ALL_WITH_H2O_LOSS+i,tmp_vals[i]-tmp_vals[half_size+i]));

	//	curr_spec_strong_inds
		
		

	//	
		const mass_t half_max = max_peak_mass*0.5;
		const mass_t half_tol = tolerance * 0.66;
		
		int num_pairs=0;
		int num_strong_pairs=0;
		float inten_both=0;
		int didx=0;
		for (i=0; i<num_peaks; i++)
		{
			const mass_t peak_mass = peaks[i].mass;
			if (peak_mass>half_max)
				break;

			const mass_t doub_mass = peak_mass * 2 - MASS_PROTON;
			const mass_t min_mass = doub_mass - half_tol;
			const mass_t max_mass = doub_mass + half_tol;
			while (didx<num_peaks && peaks[didx].mass<min_mass)
				didx++;

			if (didx==num_peaks)
				break;
			
			if (peaks[didx].mass<max_mass)
			{
				num_pairs++;
				if (curr_spec_strong_inds[i] || curr_spec_strong_inds[didx])