advancedscoremodel.cpp

MS-Clustering is designed to rapidly cluster large MS/MS datasets. The program merges similar spectr

				const int edge_idx = node.out_edge_idxs[j];
				const MultiEdge& out_edge = multi_edges[edge_idx];
				const int num_aa = out_edge.num_aa;

				if (num_aa + aa_idx >num_pep_aas)
					continue;

				const int var_idx = out_edge.get_variant_idx(num_aa,&pep_aas[aa_idx]);
				if (var_idx<0)
					continue;

				out_edge_idx = edge_idx;
				out_edge_variant = var_idx;
				break;
			}
		}


		BreakageInfo info;
		prm->fill_breakage_info(this,&info,i,in_edge_idx,in_edge_variant,out_edge_idx,out_edge_variant);
	
		node.score_combos.clear();

	//	cout << in_edge_idx << " " << in_edge_variant << " " << out_edge_idx << " " << out_edge_variant << "\t";

		info.score = score_model.score_a_single_breakage_combo(prm, node, &node.breakage, info);
	
		node.score_combos[ScoreComboLoc(info)]=info.score;
		node.score = info.score; 
		node.breakage.score = node.score;

	//	cout << node.score << endl;
	}
	prm->set_has_node_combo_scores(true);
}



/***************************************************************************
Scores the Gap-Gap combination for all nodes and sets it as the node score
****************************************************************************/
void AdvancedScoreModel::initial_combos_score(PrmGraph *prm) const
{
	const vector<int>& org_aas = config.get_org_aa();
	const vector<MultiEdge>& multi_edges = prm->get_multi_edges();
	const int num_nodes = prm->get_num_nodes();
	const mass_t mid_mass = prm->get_pm_with_19()*0.5;
	
	int i;
	for (i=0; i<num_nodes; i++)
	{
		Node& node = prm->get_non_const_node(i);
		const RegionalScoreModel& score_model = 
			regional_breakage_score_models[prm->get_charge()][prm->get_size_idx()][node.breakage.region_idx];

		vector<BreakageInfo> infos;

		BreakageInfo double_gap_info;
		prm->fill_breakage_info(this,&double_gap_info,i,NEG_INF,NEG_INF,NEG_INF,NEG_INF);
		infos.push_back(double_gap_info);

		// fill all combos for a digest node
		if (node.type == NODE_DIGEST)
		{
			int j;
			for (j=0; j<node.in_edge_idxs.size(); j++)
			{
				const int in_edge_idx = node.in_edge_idxs[j];
				const MultiEdge& in_edge = multi_edges[in_edge_idx];
				const int num_in_variants = in_edge.get_num_variants();

				int k;
				for (k=0; k<num_in_variants; k++)
				{
					int a;
					for (a=0; a<node.out_edge_idxs.size(); a++)
					{
						const int out_edge_idx = node.out_edge_idxs[a];
						const MultiEdge& out_edge = multi_edges[out_edge_idx];
						const int num_out_variants = out_edge.get_num_variants();

						int b;
						for (b=0; b<num_out_variants; b++)
						{
							BreakageInfo info;
							prm->fill_breakage_info(this,&info,i,in_edge_idx,k,out_edge_idx,b);
							if (info.connects_to_C_term || info.connects_to_N_term)
								infos.push_back(info);
						}
					}
				}
			}
		}
		
		if (node.mass>mid_mass)
		{
			int j;
			for (j=0; j<node.in_edge_idxs.size(); j++)
			{
				const int in_edge_idx = node.in_edge_idxs[j];
				const MultiEdge& in_edge = multi_edges[in_edge_idx];
				const int num_in_variants = in_edge.get_num_variants();

				int k;
				for (k=0; k<num_in_variants; k++)
				{
					BreakageInfo c_gap_info;
					prm->fill_breakage_info(this,&c_gap_info,i,in_edge_idx,k,NEG_INF,NEG_INF);
					infos.push_back(c_gap_info);
				}
			}
		}
		else
		{
			int j;
			for (j=0; j<node.out_edge_idxs.size(); j++)
			{
				const int out_edge_idx = node.out_edge_idxs[j];
				const MultiEdge& out_edge = multi_edges[out_edge_idx];
				const int num_out_variants = out_edge.get_num_variants();

				int k;
				for (k=0; k<num_out_variants; k++)
				{
					BreakageInfo n_gap_info;
					prm->fill_breakage_info(this,&n_gap_info,i,NEG_INF,NEG_INF,out_edge_idx,k);
					infos.push_back(n_gap_info);
				}
			}
		}


	

		node.score_combos.clear();

	
		int j;
		for (j=0; j<infos.size(); j++)
			infos[j].score = score_model.score_a_single_breakage_combo(prm, node, &node.breakage, infos[j]);

		for (j=0; j<infos.size(); j++)
			node.score_combos[ScoreComboLoc(infos[j])]=infos[j].score;
		
		score_t max_score=NEG_INF;
		for (j=1; j<infos.size(); j++)
			if (infos[j].score>max_score)
				max_score=infos[j].score;

		node.score = max_score; 
		node.breakage.score = node.score;
	}

	prm->set_has_node_combo_scores(true);
}

// performs scoring on demand (if the combo was not previously scored, calculates
// values, otherwise returns hashed value
score_t AdvancedScoreModel::get_node_combo_score(PrmGraph *prm, int node_idx, 
										 int in_edge_idx,  int in_var_idx, 
										 int out_edge_idx, int out_var_idx) const
{
	Config *config = prm->get_config();
	const Node& node = prm->get_node(node_idx);
	ScoreComboLoc loc(in_edge_idx,out_edge_idx,in_var_idx,out_var_idx);

	const map<ScoreComboLoc,score_t>::const_iterator it = node.score_combos.find(loc);
	if (it != node.score_combos.end())
		return it->second;

	if (node.type == NODE_N_TERM || node.type == NODE_C_TERM)
	{
		const score_t terminal_score=config->get_terminal_score();
		Node& non_const_node = prm->get_non_const_node(node_idx);
		non_const_node.score_combos[loc]=terminal_score;
		return terminal_score;
	}

	const RegionalScoreModel& score_model = 
			regional_breakage_score_models[prm->get_charge()][prm->get_size_idx()][node.breakage.region_idx];
	Node& non_const_node = prm->get_non_const_node(node_idx);

	BreakageInfo info;
	prm->fill_breakage_info(this,&info,node_idx,in_edge_idx,in_var_idx,out_edge_idx,out_var_idx);

	// score not here, calculate the combo score, store it, and return it
	const int charge = prm->get_charge();
	const int size_idx = prm->get_size_idx();
	Spectrum *spec = prm->get_source_spectrum();
	mass_t pm_with_19 = prm->get_pm_with_19();
	const RegionalScoreModel& rsm = regional_breakage_score_models[charge][size_idx][info.breakage->region_idx];

	score_t combo_score = rsm.score_a_single_breakage_combo(prm, non_const_node, info.breakage,info);

	non_const_node.score_combos[loc]=combo_score;
	non_const_node.breakage.score = combo_score;

	return combo_score;

}







/****************************************************************************
Returns the score of the constant element in the breakage:
strong features tha don't depend on aa, and all weak features
*****************************************************************************/
void RegionalScoreModel::calc_constant_element(
							   Node& node,
							   Spectrum *spec, 
							   mass_t pm_with_19,  
							   const Breakage *breakage) const
{
	node.const_strong_exps.clear();
	node.const_strong_exps.resize(strong_models.size(),0);

	int f;
	for (f=0; f<this->strong_models.size(); f++)
	{
		const StrongFragModel& strong_model = strong_models[f];
		const int frag_idx = strong_model.model_frag_idx;

		if (! strong_model.ind_has_models)
			continue;

		if (! breakage->is_frag_type_visible(frag_idx))
			continue;

		ME_Regression_Sample sam;	
		strong_model.fill_constant_vals(spec,pm_with_19,breakage,sam.f_vals);

		const int pos = breakage->get_position_of_frag_idx(frag_idx);
		if (pos>=0)
		{
			node.const_strong_exps[f]=strong_model.inten_model.get_sum_exp(sam);
		}
		else
			node.const_strong_exps[f]=strong_model.no_inten_model.get_sum_exp(sam);
	}

	node.const_regular_exps.clear();
	node.const_regular_exps.resize(regular_models.size(),0);

	for (f=0; f<regular_models.size(); f++)
	{
		const RegularFragModel& regular_model = regular_models[f];
		const int frag_idx = regular_model.model_frag_idx;

		if (! regular_model.ind_has_models)
			continue;

		if (! breakage->is_frag_type_visible(frag_idx))
			continue;

		ME_Regression_Sample sam;
		regular_model.fill_constant_vals(spec,pm_with_19,breakage,sam.f_vals);

		const int pos = breakage->get_position_of_frag_idx(frag_idx);
		if (pos>=0)
		{
			node.const_regular_exps[f] = regular_model.inten_model.get_sum_exp(sam);
		}
		else
			node.const_regular_exps[f] = regular_model.no_inten_model.get_sum_exp(sam);
	}
}


/****************************************************************************
Returns the score of the variable element in the breakages:
the strong features that are aa dependant.
*****************************************************************************/
score_t RegionalScoreModel::score_a_single_breakage_combo(
							   PrmGraph *prm,  
							   Node& node,
							   const Breakage *breakage,
							   BreakageInfo& info,
							   bool verbose) const
{
	if (node.type == NODE_N_TERM || node.type == NODE_C_TERM)
	{
		const score_t terminal_score=config->get_terminal_score();
		info.score = terminal_score;
		return terminal_score;
	}
	
	Spectrum *spec = prm->get_source_spectrum();
	const mass_t pm_with_19 = prm->get_pm_with_19();

	if (node.const_strong_exps.size()==0)
		calc_constant_element(node,spec,pm_with_19,breakage);

	score_t score=0;

	int f;
	for (f=0; f<this->strong_models.size(); f++)
	{
		const StrongFragModel& strong_model = strong_models[f];
		const int frag_idx = strong_model.model_frag_idx;

		if (! strong_model.ind_has_models)
			continue;

		if (! breakage->is_frag_type_visible(frag_idx))
			continue;

		ME_Regression_Sample sam;
		sam.f_vals.clear();
		strong_model.fill_aa_variable_vals(spec,pm_with_19,breakage,&info,sam.f_vals);
	
		score_t prev = score;
		const int pos = breakage->get_position_of_frag_idx(frag_idx);
		if (pos>=0)
		{
			const float var_exp = strong_model.inten_model.get_sum_exp(sam);
			const float e = exp(var_exp + node.const_strong_exps[f]);

			float prob = e/(1.0 + e);
			if (prob>0.99)
				prob=0.99;
			if (prob<0.001)
				prob=0.001;

	//		prob *= strong_inten_weights[f];
	//		prob += strong_inten_danc_part[f];

			const float log_random_peak = spec->get_peak(breakage->fragments[pos].peak_idx).log_rand_prob;
			score += (strong_model.inten_log_scaling_factor + log(prob) - log_random_peak);
		}
		else
		{
			const float var_exp = strong_model.no_inten_model.get_sum_exp(sam);
			const float e = exp(var_exp + node.const_strong_exps[f]);

			float prob = e/(1.0 + e);
			if (prob>0.99)
				prob=0.99;
			if (prob<0.01)
				prob=0.01;

	//		prob *= strong_no_inten_weights[f];
	//		prob += strong_no_inten_danc_part[f];


			score += (strong_model.no_inten_log_scaling_factor + log(prob) - log_one_minus_random);
		}

		if (verbose)
		{
			cout << setprecision(4) << fixed << f << "\t" << score-prev << endl;
		}
	}

	for (f=0; f<regular_models.size(); f++)
	{
		const RegularFragModel& regular_model = regular_models[f];
		const int frag_idx = regular_model.model_frag_idx;

		if (! regular_model.ind_has_models)
			continue;

		if (! breakage->is_frag_type_visible(frag_idx))
			continue;

		ME_Regression_Sample sam;
		sam.f_vals.clear();
		regular_model.fill_aa_variable_vals(spec,pm_with_19,breakage,&info,sam.f_vals);
		
		score_t prev= score;

		const int pos = breakage->get_position_of_frag_idx(frag_idx);
		if (pos>=0)
		{
			const float var_exp = regular_model.inten_model.get_sum_exp(sam);
			const float e = exp(var_exp + node.const_regular_exps[f]);

			float prob = e/(1.0 + e);
			if (prob>0.99)
				prob=0.99;
			if (prob<0.01)
				prob=0.01;

		//	prob *= regular_inten_weights[f];
		//	prob += regular_inten_danc_part[f];

			const float log_random_peak = spec->get_peak(breakage->fragments[pos].peak_idx).log_rand_prob;
			score += (regular_model.inten_log_scaling_factor + log(prob) - log_random_peak);
		}
		else
		{
			const float var_exp = regular_model.no_inten_model.get_sum_exp(sam);
			const float e = exp(var_exp + node.const_regular_exps[f]);

			float prob = e/(1.0 + e);
			if (prob>0.99)
				prob=0.99;
			if (prob<0.02)
				prob=0.02;

		//	prob *= regular_no_inten_weights[f];
		//	prob += regular_no_inten_danc_part[f];
		//	prob += 0.1;

		//	if (f>3) xxx
		//		prob = 1.0 - get_frag_prob(regular_models[f].model_frag_idx);

			score += (regular_model.no_inten_log_scaling_factor + log(prob) - log_one_minus_random);
		}

		if (verbose)
		{
			cout << f << "\t" << score-prev << endl;
		}
	}

	// correct for digest node scores
	if (node.type == NODE_DIGEST)
	{
		const score_t digest_score=config->get_digest_score();

		if ((info.connects_to_N_term && info.n_edge_idx<0) || 
			(info.connects_to_C_term && info.c_edge_idx<0) )
		{
			score -= digest_score; // penalty for ending here!
		}
			
		if (info.connects_to_N_term && info.preferred_digest_aa_N_term && score<0)
			score = 0;
		if (info.connects_to_C_term && info.preferred_digest_aa_C_term && score<0)
			score = 0;
	}

	return score;
}