pmcsqs.cpp

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

#include "PMCSQS.h"
#include "auxfun.h"


void PMCSQS_Scorer::set_sqs_mass_thresholds()
{
	sqs_mass_thresholds.clear();
	sqs_mass_thresholds.push_back(800);
	sqs_mass_thresholds.push_back(1200);
}


void PMCSQS_Scorer::set_default_sqs_correct_factors()
{
	sqs_correction_factors.resize(4);
	sqs_mult_factors.resize(4);

	sqs_correction_factors[1].push_back(0);
	sqs_correction_factors[1].push_back(-0.1);
	sqs_correction_factors[1].push_back(-0.15);

	sqs_correction_factors[2].push_back(0);
	sqs_correction_factors[2].push_back(0.05);
	sqs_correction_factors[2].push_back(0.05);

	sqs_correction_factors[2].push_back(0);
	sqs_correction_factors[2].push_back(0.05);
	sqs_correction_factors[2].push_back(0.05);

	int i;
	for (i=0; i<4; i++)
	{
		int j;
		for (j=0; j<sqs_correction_factors.size(); j++)
			sqs_mult_factors[i].push_back(1.0/(1.0+sqs_correction_factors[i][j]));
	}
}


void PMCSQS_Scorer::init_sqs_correct_factors(int max_charge, int num_sizes)
{
	sqs_correction_factors.resize(max_charge+max_charge);
	sqs_mult_factors.resize(max_charge+1);

	int i;
	for (i=1; i<=max_charge; i++)
	{
		sqs_correction_factors[i].clear();
		sqs_correction_factors[i].resize(num_sizes,0);
		sqs_mult_factors[i].clear();
		sqs_mult_factors[i].resize(num_sizes,1);
	}
}





/***************************************************************************
Gives detailed m/z and prob values for the different charges.
Returns the highest sqs probability found for a spectrum.
****************************************************************************/
float PMCSQS_Scorer::get_pmcsqs_results_for_spectrum(Config *config, 
													 const BasicSpectrum& bs,
													 vector<PmcSqsChargeRes>& res)
{
	static ME_Regression_Sample sqs_sam;

	if (! init_for_current_spec(config,bs))
		return 0; // corrupt file

	calculate_curr_spec_pmc_values(bs,bin_increment);

	const int num_charges = pmc_rank_models.size();
	const int sqs_size_idx = this->get_sqs_size_idx(bs.ssf->m_over_z);

	if (res.size()<num_charges)
		res.resize(num_charges);

	
	vector<float> min_comp_probs;		// the minimal probability when used in a comparison model
	min_comp_probs.resize(num_charges,1.0);
	min_comp_probs[0]=0;

	if (ind_initialized_sqs)
	{
		fill_fval_vector_with_SQS(bs, sqs_sam);
	
		int charge;
		for (charge=1; charge<num_charges; charge++)
		{
			
			float prob = sqs_models[charge][0][sqs_size_idx]->p_y_given_x(0,sqs_sam);

			// correct prob
			prob += sqs_correction_factors[charge][sqs_size_idx];
			prob *= sqs_mult_factors[charge][sqs_size_idx];
			if (prob<0.0)
				prob=0.0;

			res[charge].sqs_prob=prob;
		}
		
		int c;
		for (c=1; c<num_charges-1; c++)
		{
			float comp_prob = 2.0;
			int d;
			for (d=c+1; d<num_charges; d++)
			{
				if (! sqs_models[d][c][sqs_size_idx])
					continue;

				comp_prob=sqs_models[d][c][sqs_size_idx]->p_y_given_x(0,sqs_sam);

				if (comp_prob<min_comp_probs[d])
					min_comp_probs[d]=comp_prob;

				const float one_minus_prob = 1.0-comp_prob;
				if (one_minus_prob<min_comp_probs[c])
					min_comp_probs[c]=one_minus_prob;
			}
		}	
	}
	else // give the max prob charge to the input charge, the rest get 0
	{
		if (bs.ssf->charge<=0)
		{
			cout << "Error: no SQS model was read, so charge must be supplied in the spectrum!" << endl;
			exit(1);
		}

		int c;
		for (c=1; c<num_charges; c++)
			if (c != bs.ssf->charge)
				min_comp_probs[c]=0;
	}

	int charge;
	for (charge=0; charge<min_comp_probs.size(); charge++)
	{
	//	if (min_comp_probs[charge]>1.5)
	//		min_comp_probs[charge]=0;
//		cout << charge << "\t" << min_comp_probs[c] << endl;
	}


	float max_prob=0;
	for (charge=1; charge<num_charges; charge++)
	{
		

		const float prob = min_comp_probs[charge];
		if (prob>max_prob)
			max_prob=prob;

		res[charge].min_comp_prob = prob;

		if (prob>0.02)
		{
			find_best_mz_values_from_rank_model(bs,charge,res[charge]);
		}
		else // only give one guess, the second one will come from another charge
		{
			res[charge].mz1=bs.ssf->m_over_z;
			res[charge].score1=0;
			res[charge].mz2=NEG_INF;
			res[charge].score2=NEG_INF;
		}
	}
	

	return max_prob;
}


/********************************************************************
Computes the best sqs value for the spectrum. It is faster than the charge
m/z function since it only looks at a few mass positions. 
*********************************************************************
float PMCSQS_Scorer::get_sqs_for_spectrum(Config *config, 
										  const BasicSpectrum& bs, 
										  int *max_charge,
										  bool verbose)
{
	static vector<PmcSqsChargeRes> res;

	if (! this->ind_initialized_pmcr)
	{
		cout << "Error: no PMC model was read!" << endl;
		if (bs.ssf->charge <=0)
			cout << "Spectrum was read with charge <=0, so charge selection is needed!" << endl;
		exit(1);
	}

	get_pmcsqs_results_for_spectrum(config,bs,res);

	float best_prob=-1;
	int best_charge=0;

	int charge;
	for (charge =1; charge<=max_model_charge; charge++)
		if (res[charge].min_comp_prob>best_prob)
		{
			best_charge = charge;
			best_prob = res[charge].min_comp_prob;
		}

	if (max_charge)
		*max_charge = best_charge;

	return best_prob;
}*/

/********************************************************************
Computes the best sqs value for the spectrum. It is faster than the charge
m/z function since it only looks at a few mass positions. 
*********************************************************************/
float PMCSQS_Scorer::get_sqs_for_spectrum(Config *config, 
										  const BasicSpectrum& bs, 
										  int *max_charge,
										  bool verbose)
{
	static ME_Regression_Sample sqs_sam;
	static vector<QCPeak> sqs_peaks;
	static int num_sqs_peaks = 0;
	BasicSpectrum sqs_bs;

	if (verbose)
	{
		cout << "ORG spectrum: " << bs.num_peaks << endl;
//		bs.print_peaks();
		cout << endl;
	}

	if (num_sqs_peaks< 2000 || num_sqs_peaks<bs.num_peaks)
	{
		num_sqs_peaks = (int)(bs.num_peaks * 1.5);
		if (num_sqs_peaks<2000)
			num_sqs_peaks = 2000;
		
		if (num_sqs_peaks>100000)
		{
			cout << "Error: too many peak in spectrum: " << bs.num_peaks << endl;
			exit(1);
		}
		sqs_peaks.resize(num_sqs_peaks);
	}

	sqs_bs.peaks = &sqs_peaks[0];
	sqs_bs.ssf = bs.ssf;

	int best_charge=0;

	// filter peaks to acheive the required peak density for the models

	int new_num_peaks =0;
	create_filtered_peak_list_for_sqs(bs.peaks,bs.num_peaks,sqs_bs.peaks,new_num_peaks);
	sqs_bs.num_peaks= new_num_peaks;

	if (verbose)
	{
		cout << "ATER FILTERING: " << new_num_peaks << endl;
//		sqs_bs.print_peaks();
		cout << endl;
	}

	if (! init_for_current_spec(config,sqs_bs))
		return 0; // corrupt spectrum

	calculate_curr_spec_pmc_values(sqs_bs,bin_increment*3);

	const int num_charges = pmc_rank_models.size();
	const int sqs_size_idx = this->get_sqs_size_idx(sqs_bs.ssf->m_over_z);

	float max_prob=-1;
	
	fill_fval_vector_with_SQS(sqs_bs, sqs_sam);
	
	int charge;
	for (charge=1; charge<num_charges; charge++)
	{
			
		float prob = sqs_models[charge][0][sqs_size_idx]->p_y_given_x(0,sqs_sam);

		// correct prob
		prob += sqs_correction_factors[charge][sqs_size_idx];
		prob *= sqs_mult_factors[charge][sqs_size_idx];
		if (prob<0.0)
			prob=0.0;

		if (prob>max_prob)
		{
			max_prob=prob;
			best_charge=charge;
		}
	}

	if (max_charge)
		*max_charge = best_charge;

	return max_prob;
}



/******************************************************************************
Selects the the two best values of charges 1,2,3
returns the max prob found
*******************************************************************************/
float PMCSQS_Scorer::get_best_mz_charge(Config *config, const BasicSpectrum& bs, 
						   mass_t* mz1, int* charge1, float *prob1,
						   mass_t* mz2, int* charge2, float *prob2,
						   vector<PmcSqsChargeRes>* all_res)
{
	static vector<PmcSqsChargeRes> res;

	if (! this->ind_initialized_pmcr)
	{
		cout << "Error: no PMC model was read!" << endl;
		if (bs.ssf->charge <=0)
			cout << "Spectrum was read with charge <=0, so charge selection is needed!" << endl;
		exit(1);
	}

	get_pmcsqs_results_for_spectrum(config,bs,res);

	float best_prob=-1;
	int best_charge=0;

	int charge;
	for (charge =1; charge<=max_model_charge; charge++)
		if (res[charge].min_comp_prob>best_prob)
		{
			best_charge = charge;
			best_prob = res[charge].min_comp_prob;
		}

	const PmcSqsChargeRes& cr = res[best_charge];
	float second_best_prob=1.0 - (cr.score1-cr.score2)/(fabs(cr.score1)+fabs(cr.score2)+2.0);
	int   second_best_charge = best_charge;
	
	for (charge =1; charge<=max_model_charge; charge++)
		if (charge != best_charge && res[charge].min_comp_prob>second_best_prob)
		{
			second_best_charge = charge;
			second_best_prob = res[charge].min_comp_prob;
		}

	*mz1 = (mass_t)res[best_charge].mz1;
	*charge1 = best_charge;
	*prob1 = best_prob;

	if (mz2 && charge2)
	{
		*charge2 = second_best_charge;
		*prob2   = second_best_prob;
		if (second_best_charge == best_charge)
		{
			*mz2 = (mass_t)res[best_charge].mz2;
		}
		else
			*mz2 = (mass_t)res[second_best_charge].mz1;
	}

	if (all_res)
		*all_res = res;

	return best_prob;
}



/********************************************************************
Select a set of mzs nd charges that have high probabilities.
If add_one_Da_offsets is set, add -1,-2,+1,+2 mzs (if there are many
mzs, not all these offfsets will be added)
Uses emprically set probability threhsolds to chose corrected pms.
*********************************************************************/
void PMCSQS_Scorer::select_pms_and_charges(Config *config, 
								const BasicSpectrum& bs,
								vector<mass_t>& pms_with_19,
								vector<int>&    charges,
								bool	add_one_Da_offsets,
								vector<PmcSqsChargeRes>* all_res)
{
	static const mass_t offsets[]={-MASS_PROTON,MASS_PROTON};
	static const int num_all_offsets = sizeof(offsets)/sizeof(mass_t);
	static const float min_comp_prob_for_adding_second = 0.04;
	static const float min_sqs_prob_for_adding_second = 0.05;

	const int spec_charge = bs.ssf->charge;
	const mass_t spec_mz = bs.ssf->m_over_z;
	const mass_t half_tol = config->get_tolerance() * 0.5;

	pms_with_19.clear();
	charges.clear();
	int    specific_charge=0;
	mass_t specific_mz=0;
	

	if (config->get_use_spectrum_charge())
		specific_charge=spec_charge;

	static vector<PmcSqsChargeRes> res;
	int max_charge=0;
	float max_prob = 0;
	if (config->get_use_spectrum_charge() && (config->get_use_spectrum_mz() || config->get_pm_tolerance()<0.1))
	{
		pms_with_19.push_back(spec_mz * spec_charge - (spec_charge-1)*MASS_PROTON);
		charges.push_back(spec_charge);
		max_charge=spec_charge;
	}
	else
	{
		if (! this->ind_initialized_pmcr)
		{
			cout << "Error: no PMC model was read!" << endl;
			if (bs.ssf->charge <=0)
				cout << "Spectrum was read with charge <=0, so charge selection is needed!" << endl;
			exit(1);
		}

		get_pmcsqs_results_for_spectrum(config,bs,res);

		if (specific_charge>0)
		{
			max_charge = spec_charge;
			max_prob = res[spec_charge].min_comp_prob;
		}
		else
		{
			int c;
			for (c=1; c<res.size(); c++)
			{
				if (res[c].min_comp_prob>max_prob)
				{
					max_prob = res[c].min_comp_prob;
					max_charge = c;
				}
			}
		}

		pms_with_19.push_back(res[max_charge].mz1 * max_charge - (max_charge-1)*MASS_PROTON);
		charges.push_back(max_charge);
		if (res[max_charge].mz2>0)
		{
			pms_with_19.push_back(res[max_charge].mz2 * max_charge - (max_charge-1)*MASS_PROTON);
			charges.push_back(max_charge);
		}

	}

	


	// only added to the first pm_with_19, check that it doesn't overlap with others
	if (add_one_Da_offsets)
	{
		int num_offsets = num_all_offsets;
		int i;
		for (i=0; i<num_offsets; i++)
		{
			const mass_t pm_with_19 = pms_with_19[0] + offsets[i];
			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;