advancedscoremodel_regional.cpp

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

					cout << "Error: bad aa 2nd before n-side!" << endl;
					exit(1);
				}
				int aas[3]={n_second_before_cut,info->n_aa,info->c_aa};
				info->n_double_span_cat = model->get_aa_category(3,aas,info->connects_to_N_term ,info->connects_to_C_term && info->c_edge_is_single);
			}

			if (! info->c_edge_is_single)
			{
				if (c_second_after_cut<0)
				{
					cout << "Error: bad aa 2nd after c-side!" << endl;
					exit(1);
				}
				int aas[3]={info->n_aa,info->c_aa,c_second_after_cut};
				info->c_double_span_cat = model->get_aa_category(3,aas,info->connects_to_N_term && info->n_edge_is_single, info->connects_to_C_term);
			}
		}

	}
	else
	{
		info->c_aa=Gap;
		info->cc_aa=Gap;
	}
}

/************************************************************************************************
*************************************************************************************************/
void PrmGraph::extract_breakage_infos_for_score_training(Model *model,
														 int frag_idx,
														 int target_region_idx,
														 bool ind_strong_frag,
														 vector<BreakageInfo>& good_examples,
														 vector<BreakageInfo>& bad_examples) const
{

	const double Gap_ratio = 0.05;
	const mass_t tolerance = config->get_tolerance();

	const Peptide& true_pep = source_spectrum->get_peptide();
	const vector<int>& aas = true_pep.get_amino_acids();
	vector<mass_t> correct_break_masses;
	vector<int>    node_to_breakages, correct_node_idxs;
	vector<int>    correct_edge_variant_map;
	
	// find minimal and maximal nodes for which the frag is visible
	const mass_t min_mass = source_spectrum->get_min_peak_mass()-1.0;
	const mass_t max_mass = source_spectrum->get_max_peak_mass()+1.0;
	const FragmentType& frag = config->get_fragment(frag_idx);
	int min_viz_idx=nodes.size()+1;
	int max_viz_idx=0;
	int i;
	for (i=0; i<nodes.size(); i++)
	{
		mass_t exp_mass = frag.calc_expected_mass(nodes[i].mass,pm_with_19);
		if (exp_mass>min_mass && exp_mass<max_mass)
		{
			if (i<min_viz_idx)
				min_viz_idx=i;
			if (i>max_viz_idx)
				max_viz_idx=i;
		}
	}

	good_examples.clear();
	bad_examples.clear();

	true_pep.calc_expected_breakage_masses(config,correct_break_masses);
	node_to_breakages.resize(nodes.size(),NEG_INF);
	correct_node_idxs.clear();
	correct_edge_variant_map.resize(multi_edges.size(),NEG_INF);

	
	for (i=0; i<correct_break_masses.size(); i++)
	{
		const int max_node_idx = get_max_score_node(correct_break_masses[i],tolerance);
		if (max_node_idx>=0)
		{
			node_to_breakages[max_node_idx]=i;
			correct_node_idxs.push_back(max_node_idx);
		}
	}

//	cout << endl <<true_pep.as_string(config) << endl;
	for (i=0; i<multi_edges.size(); i++)
	{
		const MultiEdge& edge = multi_edges[i];
		if (node_to_breakages[edge.n_idx]>=0 && node_to_breakages[edge.c_idx]>=0)
		{
			const int brekage_idx = node_to_breakages[edge.n_idx];
			correct_edge_variant_map[i]=multi_edges[i].get_variant_idx(edge.num_aa,&aas[brekage_idx]);

			if (0 && correct_edge_variant_map[i]>=0)
			{
				cout << brekage_idx << "\t" << i << "\t" << correct_edge_variant_map[i] << "\t";
				int *var_ptr = edge.variant_ptrs[correct_edge_variant_map[i]];
				int num_aa = *var_ptr++;
				int j;
				for (j=0; j<num_aa; j++)
					cout << config->get_aa2label()[*var_ptr++];
				cout << endl;
			}
		}
	}

	vector<int> correct_in_edge_idxs;
	vector<int> correct_out_edge_idxs;
	correct_in_edge_idxs.resize(correct_node_idxs.size(),NEG_INF);
	correct_out_edge_idxs.resize(correct_node_idxs.size(),NEG_INF);

	bool had_good_connect_to_n_term = false;
	bool had_good_connect_to_c_term = false;

	// create infos for good peak samples
	for (i=0; i<correct_node_idxs.size(); i++)
	{
		const int node_idx = correct_node_idxs[i];
		const int node_region_idx = nodes[node_idx].breakage.region_idx;
		if (node_region_idx != target_region_idx || node_idx==0 || 
			node_idx==nodes.size()-1 || node_idx<min_viz_idx || node_idx>max_viz_idx)
			continue;

		const Node& node = nodes[node_idx];
		
		int n_edge_idx = NEG_INF;
		int c_edge_idx = NEG_INF;
		int n_varaint_idx = NEG_INF;
		int c_variant_idx = NEG_INF;

		int j;
		for (j=0; j<node.in_edge_idxs.size(); j++)
			if (correct_edge_variant_map[node.in_edge_idxs[j]]>=0)
			{
				n_edge_idx=node.in_edge_idxs[j];
				n_varaint_idx=correct_edge_variant_map[node.in_edge_idxs[j]];
				correct_in_edge_idxs[i]=n_edge_idx;
				break;
			}

		for (j=0; j<node.out_edge_idxs.size(); j++)
			if (correct_edge_variant_map[node.out_edge_idxs[j]]>=0)
			{
				c_edge_idx=node.out_edge_idxs[j];
				c_variant_idx=correct_edge_variant_map[node.out_edge_idxs[j]];
				correct_out_edge_idxs[i]=c_edge_idx;
				break;
			}

		BreakageInfo info;
		fill_breakage_info(model,&info,node_idx,n_edge_idx,n_varaint_idx,c_edge_idx,c_variant_idx,1);

		if (info.connects_to_N_term)
			had_good_connect_to_n_term=true;
		
		if (info.connects_to_C_term)
			had_good_connect_to_c_term=true;

		good_examples.push_back(info);
		if (n_edge_idx>=0 && c_edge_idx>=0 && my_random()<Gap_ratio*0.33)
		{
			BreakageInfo gap_info;
			fill_breakage_info(model,&gap_info,node_idx,NEG_INF,NEG_INF,c_edge_idx,c_variant_idx,11);
			good_examples.push_back(gap_info);
		}

		if (n_edge_idx>=0 && c_edge_idx>=0 && my_random()<Gap_ratio*0.33)
		{
			BreakageInfo gap_info;
			fill_breakage_info(model,&gap_info,node_idx,n_edge_idx,n_varaint_idx,NEG_INF,NEG_INF,11);
			good_examples.push_back(gap_info);
		}
	}

	// create for bad peak samples
	const int num_good = good_examples.size();

	// same as good nodes, but using bad edges
	// perform only for strong fragments!
	const int num_half_bad_examples = (ind_strong_frag ? int(num_good*0.333) : 0);
	for (i=0; i<num_half_bad_examples; i++)
	{
		const int node_idx = correct_node_idxs[i];
		const int node_region_idx = nodes[node_idx].breakage.region_idx;

		if (node_region_idx != target_region_idx || node_idx==0 || 
			node_idx==nodes.size()-1 || node_idx<min_viz_idx || node_idx>max_viz_idx)
			continue;

		const Node& node = nodes[node_idx];
		if (node.breakage.get_position_of_frag_idx(frag_idx)<0) // don't use such samples if there is no peak
			continue;

		if (node.in_edge_idxs.size()>1 && correct_in_edge_idxs[i]>=0 &&
			node.out_edge_idxs.size()>1 && correct_out_edge_idxs[i]>=0)
		{
			int j;
			vector<int> in_idxs,out_idxs;
			for (j=0; j<node.in_edge_idxs.size(); j++)
			{
				const int edge_idx = node.in_edge_idxs[j];
				if (correct_edge_variant_map[edge_idx]>=0)
					continue;
				if (multi_edges[edge_idx].num_aa == 1)
				{
					in_idxs.push_back(edge_idx);
				}
				else if (multi_edges[edge_idx].num_aa != 1 && (in_idxs.size()<3 || my_random()<0.1))
					in_idxs.push_back(edge_idx);
			}
			
			for (j=0; j<node.out_edge_idxs.size(); j++)
			{
				const int edge_idx = node.out_edge_idxs[j];
				if (correct_edge_variant_map[edge_idx]>=0)
					continue;
				if (multi_edges[edge_idx].num_aa == 1)
				{
					out_idxs.push_back(edge_idx);
				}
				else if (multi_edges[edge_idx].num_aa != 1 && (out_idxs.size()<3 || my_random()<0.1))
					out_idxs.push_back(edge_idx);
			}

			if (in_idxs.size()==0 || out_idxs.size()==0)
				continue;

			const int bad_n_edge = in_idxs[int(my_random()*in_idxs.size())];
			const int bad_c_edge = out_idxs[int(my_random()*out_idxs.size())];
			const int n_var_idx = int(multi_edges[bad_n_edge].variant_ptrs.size() * my_random());
			const int c_var_idx = int(multi_edges[bad_c_edge].variant_ptrs.size() * my_random());

			BreakageInfo info;
			fill_breakage_info(model,&info,node_idx,bad_n_edge,n_var_idx,bad_c_edge,c_var_idx,2);
			bad_examples.push_back(info);
			if (bad_n_edge>=0 && bad_c_edge>=0 && my_random()<Gap_ratio)
			{
				BreakageInfo gap_info;
				fill_breakage_info(model,&gap_info,node_idx,NEG_INF,NEG_INF,bad_c_edge,c_var_idx,22);
				bad_examples.push_back(gap_info);
			}

			if (bad_n_edge>=0 && bad_c_edge>=0 && my_random()<Gap_ratio)
			{
				BreakageInfo gap_info;
				fill_breakage_info(model,&gap_info,node_idx,bad_n_edge,n_var_idx,NEG_INF,NEG_INF,22);
				bad_examples.push_back(gap_info);
			}
		}
	}


	// mostly single edges from the highest scoring incorrect nodes
	// gives advantage for aa combinations with high scoring nodes
	vector<score_pair> pairs;
	for (i=1; i<nodes.size()-1; i++)
		if (node_to_breakages[i]<0 && i>=min_viz_idx && i<=max_viz_idx &&
			nodes[i].breakage.region_idx==target_region_idx )
				pairs.push_back(score_pair(i,nodes[i].score));
	
	sort(pairs.begin(),pairs.end());

	vector<int> selected_idxs, random_idxs;
	int num_to_select = int(num_good * 3);
	for (i=0; i<num_to_select && i<pairs.size(); i++)
		selected_idxs.push_back(pairs[i].idx);

	if (pairs.size()< num_to_select)
		num_to_select = pairs.size();