advancedscoremodel_regional.cpp

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

			cout << "Problem with models!!!" << endl;
			exit(1);
		}

		cout << endl << endl << "TRAINING NO INTENSITY MODEL FOR CHARGE " <<
			charge << " SIZE " << size_idx << " REGION " << region_idx << " FRAGMENT " << i << " " <<
			config->get_fragment(frag_idx).label << endl << endl;

		no_inten_ds.purge_low_count_features(min_num_of_samples_per_feature);
		int num_bad=no_inten_ds.check_samples(true);
			if (num_bad>0)
				cout << "Warning: had " << num_bad << " bad samples removed!" << endl;

		no_inten_ds.print_summary();
		no_inten_ds.print_feature_summary(cout, ScoreModelFields_SNI_names);
		if (! strong_models[i].no_inten_model.train_cg(no_inten_ds,num_me_rounds,2E-5))
		{
			cout << "Coudln't train no inten ME model, exiting!" << endl;
			exit(1);
		}
		else
		{
			cout << endl << "NO INTENSITY - Charge " << charge << " size " << size_idx << " region " << region_idx << 
				" fragment " << config->get_fragment(frag_idx).label << endl;
			strong_models[i].no_inten_model.print_ds_probs(no_inten_ds);
			strong_models[i].no_inten_log_scaling_factor = 
				log(strong_models[i].no_inten_model.calc_log_scaling_constant(0,no_inten_ds,1.1*(1.0-frag_prob)));
		}

		strong_models[i].ind_has_models = true;
	}

	for (i=0; i<regular_models.size(); i++)
	{
		ME_Regression_DataSet inten_ds, no_inten_ds;
		const int frag_idx = regular_models[i].model_frag_idx;
		const float frag_prob = get_frag_prob(frag_idx);

		cout << endl << endl << "TRAINING INTENSITY MODEL FOR CHARGE " <<
			charge << " SIZE " << size_idx << " REGION " << region_idx << " FRAGMENT " << i << " " <<
			config->get_fragment(frag_idx).label << endl << endl;

		int j;
		for (j=0; j<3; j++)
		{
			cout << "SEED: " << get_random_seed() << endl;
			create_training_set(model, regular_models[i], fm, inten_ds, no_inten_ds);

			inten_ds.purge_low_count_features(min_num_of_samples_per_feature);
			int num_bad=inten_ds.check_samples(true);
			if (num_bad>0)
				cout << "Warning: had " << num_bad << " bad samples removed!" << endl;


			inten_ds.print_summary();
			inten_ds.print_feature_summary(cout, ScoreModelFields_RI_names);
			if (! regular_models[i].inten_model.train_cg(inten_ds,num_me_rounds,2E-5))
			{
				cout << "Coudln't train ME model, setting all weights to 0! (" <<j<<")"<< endl;
			}
			else
			{
				cout << endl << "INTENSTY - Charge " << charge << " size " << size_idx << " region " << region_idx << 
					" fragment " << config->get_fragment(frag_idx).label << endl ;
				regular_models[i].inten_model.print_ds_probs(inten_ds);
				regular_models[i].inten_log_scaling_factor = 
					log(regular_models[i].inten_model.calc_log_scaling_constant(0,inten_ds,1.1*frag_prob));
				break;
			}
		}
		
		cout << endl << endl << "TRAINING NO INTENSITY MODEL FOR CHARGE " <<
				charge << " SIZE " << size_idx << " REGION " << region_idx << " FRAGMENT " << i << " " <<
				config->get_fragment(frag_idx).label << endl << endl;

		no_inten_ds.purge_low_count_features(min_num_of_samples_per_feature);

		int num_bad=no_inten_ds.check_samples(true);
		if (num_bad>0)
			cout << "Warning: had " << num_bad << " bad samples removed!" << endl;


		no_inten_ds.print_summary();
		no_inten_ds.print_feature_summary(cout, ScoreModelFields_RNI_names);
		if (! regular_models[i].no_inten_model.train_cg(no_inten_ds,num_me_rounds,2E-5))
		{
			cout << "Coudln't train ME model, setting all weights to 0!" << endl;
		}
		else
		{
			cout << endl << "NO INTENSTY - Charge " << charge << " size " << size_idx << " region " << region_idx << 
				" fragment " << config->get_fragment(frag_idx).label << endl;
				regular_models[i].no_inten_model.print_ds_probs(no_inten_ds);
				regular_models[i].no_inten_log_scaling_factor = 
					log(regular_models[i].no_inten_model.calc_log_scaling_constant(0,no_inten_ds,(1.0-frag_prob)));
		}

		regular_models[i].ind_has_models = true;
	}	


	was_initialized = true;

	write_regional_score_model(name);

	return true;
}




void BreakageInfo::print(Config *config) const
{
	const vector<string>& aa2label = config->get_aa2label();
	cout << setw(2) << this->type << "> ";
	cout << (this->connects_to_N_term ? '[' : ' ');
	cout << (this->n_edge_is_single ? '1' : '2');
	cout <<  " " << aa2label[this->n_aa] << " : " << aa2label[this->c_aa] << " ";
	cout << (this->c_edge_is_single ? '1' : '2');
	cout << (this->connects_to_C_term ? ']' : ' ');
	cout << " " << fixed << setprecision(3) << "  # " <<this->node_idx << " ";
	if (this->missed_cleavage)
		cout << "MC!";
	cout << endl;
}

void PrmGraph::fill_breakage_info(const Model *model, BreakageInfo *info, int node_idx,
								  int n_edge_idx, int n_variant_idx,
								  int c_edge_idx, int c_variant_idx, int type) const
{
	const vector<mass_t>& aa2mass = config->get_aa2mass();
	const vector<int>& n_term_digest_aas = config->get_n_term_digest_aas();
	const vector<int>& c_term_digest_aas = config->get_c_term_digest_aas();
	
	info->node_idx = node_idx;
	info->breakage = &nodes[node_idx].breakage;
	info->type = type;

	info->n_edge_idx = n_edge_idx;
	info->n_var_idx  = n_variant_idx;
	info->c_edge_idx = c_edge_idx;
	info->c_var_idx  = c_variant_idx;

	int n_second_before_cut = NEG_INF;
	if (n_edge_idx>=0)
	{
		const MultiEdge& n_edge = multi_edges[n_edge_idx];
		const Node& n_node = nodes[n_edge.n_idx];
		int *v_ptr = n_edge.variant_ptrs[n_variant_idx];

		info->n_var_ptr = v_ptr;

		const int num_aa = *v_ptr++;
		int n_aa = v_ptr[num_aa-1];
		mass_t exp_edge_mass=0;
		int j;
		for (j=0; j<num_aa; j++)
			exp_edge_mass+=aa2mass[v_ptr[j]];

		if (n_aa == Ile || n_aa == Xle)
			n_aa = Leu;

		info->n_aa = n_aa;
		info->nn_aa = v_ptr[0];
		
		if (info->nn_aa == Ile || info->nn_aa == Xle)
			info->nn_aa = Leu;

		info->n_break = n_edge.n_break;
		info->exp_n_edge_mass = exp_edge_mass;
		info->n_edge_is_single = (num_aa ==1);

		if (! info->n_edge_is_single)
			n_second_before_cut = v_ptr[num_aa-2];
		if (n_second_before_cut == Ile || n_second_before_cut == Xle)
			n_second_before_cut = Leu;

		if (n_node.type == NODE_N_TERM)
		{
			info->connects_to_N_term=true;
			int j;
			for (j=0; j<n_term_digest_aas.size(); j++)
				if (n_term_digest_aas[j] == *v_ptr)
				{
					info->preferred_digest_aa_N_term=true;
					break;
				}
		}

		for (j=0; j<c_term_digest_aas.size(); j++)
			if (c_term_digest_aas[j] == v_ptr[num_aa-1])
			{
				info->missed_cleavage= true;
				break;
			}

		info->ind_n_edge_overlaps = (n_edge.ind_edge_overlaps);
		info->n_side_cat = model->get_aa_category(num_aa,v_ptr,info->connects_to_N_term,false);
	}
	else
	{
		info->n_aa=Gap;
		info->nn_aa=Gap;
	}

	if (c_edge_idx>=0)
	{
		const MultiEdge& c_edge = multi_edges[c_edge_idx];
		const Node& c_node = nodes[c_edge.c_idx];
		int *v_ptr = c_edge.variant_ptrs[c_variant_idx];

		info->c_var_ptr = v_ptr;

		const int num_aa = *v_ptr++;
		int c_aa = v_ptr[0];
		mass_t exp_edge_mass=0;
		int j;
		for (j=0; j<num_aa; j++)
			exp_edge_mass+=aa2mass[v_ptr[j]];

		if (c_aa == Ile || c_aa == Xle)
			c_aa = Leu;

		info->c_aa = c_aa;
		info->cc_aa = v_ptr[num_aa-1];

		if (info->cc_aa == Ile || info->cc_aa == Xle)
			info->cc_aa = Leu;

		info->c_break = c_edge.c_break;
		info->exp_c_edge_mass = exp_edge_mass;
		info->c_edge_is_single = (num_aa ==1);

		int c_second_after_cut=NEG_INF;
		if (! info->c_edge_is_single)
			c_second_after_cut = v_ptr[1];
		if (c_second_after_cut == Ile || c_second_after_cut == Xle)
			c_second_after_cut = Leu;

		if (c_node.type == NODE_C_TERM)
		{
			info->connects_to_C_term=true;
			int j;
			for (j=0; j<c_term_digest_aas.size(); j++)
				if (c_term_digest_aas[j] == v_ptr[num_aa-1])
				{
					info->preferred_digest_aa_C_term=true;
					break;
				}
		}

		for (j=0; j<n_term_digest_aas.size(); j++)
			if (n_term_digest_aas[j] == v_ptr[0])
			{
				info->missed_cleavage= true;
				break;
			}

		info->ind_c_edge_overlaps = (c_edge.ind_edge_overlaps);
		info->c_side_cat = model->get_aa_category(num_aa,v_ptr,false,info->connects_to_C_term);
		if (n_edge_idx>=0)
		{
			int aas[2]={info->n_aa,info->c_aa};
			info->span_cat = model->get_aa_category(2,aas,info->connects_to_N_term && info->n_edge_is_single,info->connects_to_C_term && info->c_edge_is_single);
			
			if (! info->n_edge_is_single)
			{
				if (n_second_before_cut<0)
				{