advancedscoremodel_regional.cpp

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

		regular_no_inten_weights[0]=0.8;
	}

	if (num_regular>3)
	{
		regular_inten_weights[0]=0.9;
		regular_no_inten_weights[0]=0.8;
	}

	if (num_regular>4)
	{
		regular_inten_weights[0]=0.9;
		regular_no_inten_weights[0]=0.8;
	}*/

	for (i=0; i<num_regular; i++)
	{
		const float frag_prob = get_frag_prob(regular_models[i].model_frag_idx);
		regular_inten_danc_part[i]=(1.0-regular_inten_weights[i])*frag_prob;
		regular_no_inten_danc_part[i]=(1.0-regular_no_inten_weights[i])*(1.0-frag_prob);

	//	if (size_idx == 1 && region_idx == 1)
	//		cout << i << "\t" << frag_prob << "\t" << regular_inten_danc_part[i] << "\t" << regular_no_inten_danc_part[i] << endl;
	}
}



/***************************************************************************
****************************************************************************/
void RegionalScoreModel::create_training_set(Model *model,
											 const FragModel& frag_model,
											 const FileManager& fm,
											 ME_Regression_DataSet& inten_ds,
											 ME_Regression_DataSet& no_inten_ds) const
{
	const mass_t min_mass = config->get_min_mass_for_size_idx(charge,size_idx);
	const mass_t max_mass = config->get_max_mass_for_size_idx(charge,size_idx);
	const mass_t tolerance_diff = config->get_tolerance()*0.5;
	const int frag_idx = frag_model.model_frag_idx;

	bool is_strong=false;
	int s;
	for (s=0; s<strong_models.size(); s++)
		if (frag_idx == strong_models[s].model_frag_idx)
			is_strong=true;

	FileSet fs;
	fs.select_files_in_mz_range(fm,min_mass/charge-2.0,max_mass/charge+2.0,charge);
	const vector<SingleSpectrumFile *>& all_ssfs = fs.get_ssf_pointers();

	cout << "Selected " << all_ssfs.size() << " headers..." << endl;
	cout << "Min m/z " << min_mass/charge << endl;
	cout << "Max m/z " << max_mass/charge << endl;

	PMCSQS_Scorer *pmcsqs = (PMCSQS_Scorer *)model->get_pmcsqs_ptr();
	const bool use_pmcr = (pmcsqs && pmcsqs->get_ind_initialized_pmcr());

	if (use_pmcr)
		cout << "Using parent mass correction model to set PM.." << endl;

	BasicSpecReader bsr;
	vector<QCPeak> peaks;
	peaks.resize(10000);
	int i;
	for (i=0; i<all_ssfs.size(); i++)
	{
		BasicSpectrum bs;
		Spectrum s;
		PrmGraph prm;

		const mass_t true_mass_with_19 = all_ssfs[i]->peptide.get_mass_with_19();
		const int num_peaks = bsr.read_basic_spec(config,fm,all_ssfs[i],&peaks[0]);
		bs.num_peaks = num_peaks;
		bs.ssf = all_ssfs[i];
		bs.peaks = &peaks[0];

		// calc corrected pm_with_19, if it is good use it, otherwise, use a value with +- U[0,toleance/2]
		mass_t pm_with_19=NEG_INF;
	
		if (use_pmcr)
		{
			mass_t	mz1,mz2;
			int		charge1,charge2;
			float	prob1,prob2;
					

			// output m/z and prob values for the different charge states
			model->get_best_mz_charge(config,bs,&mz1,&charge1,&prob1,&mz2,&charge2,&prob2);

			const mass_t corr1_pm_with_19 = mz1*charge1 - MASS_PROTON*(charge1-1);
			const mass_t corr2_pm_with_19 = mz2*charge2 - MASS_PROTON*(charge2-1);
			if (fabs(corr2_pm_with_19-true_mass_with_19)<tolerance_diff)
				pm_with_19 = corr2_pm_with_19;
			if (fabs(corr1_pm_with_19-true_mass_with_19)<tolerance_diff)
				pm_with_19 = corr1_pm_with_19;
		}

		if (pm_with_19<0) // use a random value
		{
			double r=my_random();
			mass_t offset = r*r*tolerance_diff;
			if (my_random()<0.5)
				offset *= -1;
			pm_with_19 = true_mass_with_19 + offset;
	
		}

		int spec_size_idx = config->calc_size_idx(charge,pm_with_19);
		if (spec_size_idx != size_idx)
			continue;
	
		s.init_from_QCPeaks(config,&peaks[0],num_peaks,all_ssfs[i]);
		s.set_corrected_pm_with_19(pm_with_19);
		prm.create_graph_from_spectrum(model,&s,pm_with_19,all_ssfs[i]->charge);

	//	cout << i << endl;

	//	s.print_expected_by();
	//	s.print_spectrum();
	//	prm.print_with_combo_tables();
	//	prm.print();
	//	exit(0);

		vector<BreakageInfo> good_peak_examples, bad_peak_examples;
		prm.extract_breakage_infos_for_score_training(model,
													  frag_model.model_frag_idx, 
													  region_idx,
													  is_strong,
													  good_peak_examples, 
													  bad_peak_examples);

	//	cout << "GOOD:" << endl;
		int j;
		for (j=0; j<good_peak_examples.size(); j++)
		{
		//	good_peak_examples[j].print(config);
			ME_Regression_Sample sam;
			sam.label=0;

			frag_model.fill_single_frag_vector(&s, pm_with_19, good_peak_examples[j].breakage,
				good_peak_examples[j], sam.f_vals);

			if (good_peak_examples[j].breakage->get_position_of_frag_idx(frag_idx)>=0)
			{
				inten_ds.add_sample(sam);
			//	sam.print((is_strong ? ScoreModelFields_SI_names : ScoreModelFields_RI_names));
			}
			else
				no_inten_ds.add_sample(sam);
		}
		
	//	cout << "BAD:" << endl;
		for (j=0; j<bad_peak_examples.size(); j++)
		{
		//	bad_peak_examples[j].print(config);
			ME_Regression_Sample sam;
			sam.label=1;

			frag_model.fill_single_frag_vector(&s, pm_with_19, bad_peak_examples[j].breakage,
				bad_peak_examples[j], sam.f_vals);

			if (bad_peak_examples[j].breakage->get_position_of_frag_idx(frag_idx)>=0)
			{
				inten_ds.add_sample(sam);
			}
			else
				no_inten_ds.add_sample(sam);

	//		sam.print();
		}

		if (i>0 && i %1000 == 0)
		{
			cout << i << "/" << all_ssfs.size() << " ..." << endl;
		}
	}

	
	const score_t frag_prob  = this->get_frag_prob(frag_idx);
	cout << "Probability of observeing fragment: " << frag_prob << endl;

	inten_ds.num_classes=2;
	inten_ds.num_features = (is_strong ? (int)SI_NUM_FEATURES : (int)RI_NUM_FEATURES);
	inten_ds.calibrate_class_weights((is_strong ? 0.45 : 0.2));
	
	no_inten_ds.num_classes=2;
	no_inten_ds.num_features = (is_strong ? (int)SNI_NUM_FEATURES : (int)RNI_NUM_FEATURES);
	no_inten_ds.calibrate_class_weights((is_strong ? 0.2 : 0.5));

	if (is_strong)
	{
		vector<int> inten_features;
		inten_features.push_back(SI_IND_CONNECTS_TO_N_TERM);
		inten_features.push_back(SI_IND_N_IS_GAP);
		inten_features.push_back(SI_IND_C_IS_GAP);

		inten_ds.serial_scale(inten_features);

		vector<int> no_inten_features;
		no_inten_features.push_back(SNI_IND_CONNECTS_TO_N_TERM);
		no_inten_features.push_back(SNI_IND_N_IS_GAP);
		no_inten_features.push_back(SNI_IND_C_IS_GAP);

		no_inten_ds.serial_scale(no_inten_features);
	}
	else
	{
		vector<int> inten_features;
		inten_features.push_back(RI_IND_N_IS_GAP);
		inten_features.push_back(RI_IND_C_IS_GAP);

		inten_ds.serial_scale(inten_features);

		vector<int> no_inten_features;
		no_inten_features.push_back(RNI_IND_N_IS_GAP);
		no_inten_features.push_back(RNI_IND_C_IS_GAP);

		no_inten_ds.serial_scale(no_inten_features);
	}
}


bool RegionalScoreModel::train_regional_score_model(Model *model, const char *name, const FileManager& fm)
{
	const int min_num_of_samples_per_feature=6;
	const int num_me_rounds = 750;
	int i;
	for (i=0; i<strong_models.size(); i++)
	{
		ME_Regression_DataSet inten_ds, no_inten_ds;
		const int frag_idx = strong_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, strong_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_SI_names);
			if (! strong_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 " << frag_idx << " " << config->get_fragment(frag_idx).label <<  endl;
				strong_models[i].inten_model.print_ds_probs(inten_ds);
				strong_models[i].inten_log_scaling_factor = 
					log(strong_models[i].inten_model.calc_log_scaling_constant(0,inten_ds,1.2*frag_prob));

				break;
			}
		}

		if (j==3)
		{