peakrankpartialdenovo.cpp

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

			real_feature_names.push_back(label);
		}
	}
}



void PeakRankModel::set_partial_denovo_feature_names()
{
	const int num_aas = this->get_num_model_aas();
	int i;

	real_feature_names.clear();
	real_feature_stage_idxs.clear();

	real_feature_names.push_back("N RKH LEVEL");
	real_feature_names.push_back("C RKH LEVEL");

	real_feature_names.push_back("PROP PEAK MASS IN VIS RANGE"); 
	push_back_all_RHK_pairs(real_feature_names,"PROP PEAK MASS IN VIS RANGE");

	real_feature_names.push_back("DIS FROM MIN");
	
	real_feature_names.push_back("DIS FROM MAX");

	real_feature_names.push_back("CUT POS");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back("CUT -1 = " + model_aa_labels[i] + ", CUT POS");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back("CUT -2 = " + model_aa_labels[i] + ", CUT POS");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back("CUT -3 = " + model_aa_labels[i] + ", CUT POS");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back("CUT +1 = " + model_aa_labels[i] + ", CUT POS");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back("CUT +2 = " + model_aa_labels[i] + ", CUT POS");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back("CUT +3 = " + model_aa_labels[i] + ", CUT POS");

	for (i=0; i<num_aas; i++)
	{
		int j;
		for (j=0; j<num_aas; j++)
			real_feature_names.push_back("CUT IS " + model_aa_labels[i] + "-" + model_aa_labels[j] + ", CUT POS");
	}


	for (i=0; i<num_aas; i++)
		real_feature_names.push_back("CUT -1 = " + model_aa_labels[i] + ", RHK LEVEL");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back("CUT -2 = " + model_aa_labels[i] + ", RHK LEVEL");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back("CUT -3 = " + model_aa_labels[i] + ", RHK LEVEL");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back("CUT +1 = " + model_aa_labels[i] + ", RHK LEVEL");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back("CUT +2 = " + model_aa_labels[i] + ", RHK LEVEL");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back("CUT +3 = " + model_aa_labels[i] + ", RHK LEVEL");


	for (i=0; i<num_aas; i++)
	{
		int j;
		for (j=0; j<num_aas; j++)
			real_feature_names.push_back(model_aa_labels[i] + model_aa_labels[j] + " BEFORE CUT, CUT POS");
	}

	for (i=0; i<num_aas; i++)
	{
		int j;
		for (j=0; j<num_aas; j++)
			real_feature_names.push_back(model_aa_labels[i] + model_aa_labels[j] + " AFTER CUT, CUT POS");
	}


	for (i=0; i<num_aas; i++)
		real_feature_names.push_back( "|"+ model_aa_labels[i] + "XK, CUT POS");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back( "|"+ model_aa_labels[i] + "XXK, CUT POS");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back( "|"+ model_aa_labels[i] + "XXXK, CUT POS");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back( "|"+ model_aa_labels[i] + "XXXXK, CUT POS");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back( model_aa_labels[i] + "|XK, CUT POS");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back( model_aa_labels[i] + "|XXK, CUT POS");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back( model_aa_labels[i] + "|XXXK, CUT POS");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back( model_aa_labels[i] + "|XXXXK, CUT POS");



	for (i=0; i<num_aas; i++)
		real_feature_names.push_back( "|"+ model_aa_labels[i] + "XR, CUT POS");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back( "|"+ model_aa_labels[i] + "XXR, CUT POS");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back( "|"+ model_aa_labels[i] + "XXXR, CUT POS");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back( "|"+ model_aa_labels[i] + "XXXXR, CUT POS");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back( model_aa_labels[i] + "|XR, CUT POS");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back( model_aa_labels[i] + "|XXR, CUT POS");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back( model_aa_labels[i] + "|XXXR, CUT POS");

	for (i=0; i<num_aas; i++)
		real_feature_names.push_back( model_aa_labels[i] + "|XXXXR, CUT POS");
}




/****************************************************************************
Reduces a TrainingPeptide in size to form a denovo training peptide.
The shorter peptide might have intensities in cuts 0, and cut |length|.
*****************************************************************************/
void convert_tps_to_partial_denovo(Config *config, 
								   vector<TrainingPeptide>& all_tps, 
								   int num_to_add,
								   bool verbose)
{
	const int y_frag_idx = config->get_frag_idx_from_label("y");
	const int b_frag_idx = config->get_frag_idx_from_label("b");
	const int y2_frag_idx = config->get_frag_idx_from_label("s2+10.2");
	const int b2_frag_idx = config->get_frag_idx_from_label("b2");
	const vector<mass_t>& aa2mass = config->get_aa2mass();
	const int org_num_tps = all_tps.size();
	vector<double> rand_threshes;
	
	rand_threshes.resize(100,0);
	int i;
	for (i=0; i<=6; i++)
		rand_threshes[i]=1.0;
	rand_threshes[7] = 0.7;    rand_threshes[8] = 0.65;  rand_threshes[9]=0.6;
	rand_threshes[10]= 0.6;    rand_threshes[11]= 0.55; rand_threshes[12]=0.5;
	rand_threshes[13]= 0.45;   rand_threshes[14]=0.4;   rand_threshes[15]=0.3;
	rand_threshes[16]= 0.2;    rand_threshes[17]=0.15;  rand_threshes[18]=0.1;
	rand_threshes[19]= 0.075;  rand_threshes[20]=0.05;

	if (num_to_add>0)
		all_tps.resize(org_num_tps*(1+num_to_add));

	int tp_idx;
	for (tp_idx=0; tp_idx<org_num_tps; tp_idx++)
	{
		const TrainingPeptide org_tp = all_tps[tp_idx];
		const int org_length = org_tp.length;

		if (verbose)
		{
			cout << tp_idx << "\t";
			org_tp.print(config);
			cout << endl;
		}

		int i;
		for (i=0; i<= num_to_add; i++)
		{
			TrainingPeptide& new_tp = all_tps[i*org_num_tps + tp_idx];

			if (rand_threshes[org_length]>0 && my_random()<rand_threshes[org_length])
			{
				new_tp = org_tp;
				continue;
			}

			// create a ruff score vector to use for subset selection
			vector<int> ruff_scores;
			ruff_scores.resize(org_tp.length+1,0);
			ruff_scores[0]=2;
			ruff_scores[org_tp.length]=1;
			
			const int last_aa = org_tp.amino_acids[org_tp.length-1];
			if (last_aa==Arg || last_aa==Lys)
				ruff_scores[org_tp.length]=2;
			
			vector<int> frags;
			frags.push_back(y_frag_idx);
			frags.push_back(b_frag_idx);
			if (org_tp.charge == 3 || (org_tp.charge == 2 && org_tp.pm_with_19 > 1800))
			{
				frags.push_back(y2_frag_idx);
				frags.push_back(b2_frag_idx);
			}

			int f,c;
			for (f=0; f<frags.size(); f++)
			{
				const int pos = org_tp.get_frag_idx_pos(frags[f]);
				if (pos<0)
					continue;

				int c;
				for (c=1; c<org_tp.length; c++)
					if (org_tp.intens[pos][c]>0)
						ruff_scores[c]++;
			}
			for (c=1; c<ruff_scores.size(); c++)
				ruff_scores[c]+=ruff_scores[c-1];


			int min_length = 6;
			int max_length = org_tp.length - 3;
			if (org_tp.length == 7)
				max_length = 6;
			if (org_tp.length == 8)
				max_length = 7;
			if (org_tp.length == 9)
				max_length = 7;

			if (max_length>20)
				max_length=20;

			const int selected_length = min_length + (int)(my_random()*(max_length-min_length+1));
			
			// find best subset of this length
			const int max_start = org_tp.length - selected_length;
			int best_idx=-1;
			int best_ruff_score=0;
			

			for (c=0; c<=max_start; c++)
			{
				int score = ruff_scores[c+selected_length];
				if (c>0)
					score -= ruff_scores[c-1];
				if (score>=best_ruff_score)
				{
					best_idx=c;
					best_ruff_score = score;
				}
				if (verbose)
				{
					cout << c << ":" << score << " ";
				}
			}
			if (verbose)
				cout << endl;

			if (best_idx<0)
			{
				cout << "Error: something is wrong with the ruff scores!" << endl;
				exit(1);
			}

			// create new tp
			const vector<int>& amino_acids = org_tp.amino_acids;
			int j;
			new_tp.amino_acids.resize(selected_length,0);
			mass_t n_mass =0;
			for (j=0; j<best_idx; j++)
				n_mass+=aa2mass[amino_acids[j]];

			new_tp.n_mass = n_mass;
			mass_t total_mass = n_mass;
			for (j=0; j<selected_length; j++)
			{
				new_tp.amino_acids[j]=amino_acids[j+best_idx];
				total_mass += aa2mass[new_tp.amino_acids[j]];
			}

			if (total_mass>org_tp.pm_with_19)
			{
				cout << "Error: mismatch in masses with de novo tp!" << endl;
				exit(1);
			}



			if (best_idx>0)
			{
				int j;
				for (j=0; j<best_idx; j++)
					if (amino_acids[j] == His)
						new_tp.best_n_removed = His;
				for (j=0; j<best_idx; j++)
					if (amino_acids[j] == Lys)
						new_tp.best_n_removed = Lys;
				for (j=0; j<best_idx; j++)
					if (amino_acids[j] == Arg)
						new_tp.best_n_removed = Arg;
			}

			const int c_idx = best_idx + selected_length;
			if (c_idx<org_tp.length)
			{
				int j;
				for (j=c_idx; j<org_tp.length; j++)
					if (amino_acids[j] == His)
						new_tp.best_c_removed = His;
				for (j=c_idx; j<org_tp.length; j++)
					if (amino_acids[j] == Lys)
						new_tp.best_c_removed = Lys;
				for (j=c_idx; j<org_tp.length; j++)
					if (amino_acids[j] == Arg)
						new_tp.best_c_removed = Arg;
			}

			new_tp.length = selected_length;
			new_tp.charge = org_tp.charge;
			new_tp.frag_idxs = org_tp.frag_idxs;
			
			// change mobility according to observed + removed aas
			int num_arg=0,num_lys=0,num_his=0;
			for (j=0; j<amino_acids.size(); j++)
			{
				if (amino_acids[j]==Arg)
					num_arg++;
				if (amino_acids[j]==Lys)
					num_lys++;
				if (amino_acids[j]==His)
					num_his++;
			}

			if (new_tp.best_c_removed>0)
			{
				if (new_tp.best_c_removed == Arg) num_arg++;
				if (new_tp.best_c_removed == Lys) num_lys++;
				if (new_tp.best_c_removed == His) num_his++;
			}

			if (new_tp.best_n_removed>0)
			{
				if (new_tp.best_n_removed == Arg) num_arg++;
				if (new_tp.best_n_removed == Lys) num_lys++;
				if (new_tp.best_n_removed == His) num_his++;
			}

			new_tp.mobility = get_proton_mobility(org_tp.charge,num_arg,num_his,num_lys);
			
			new_tp.pm_with_19 = org_tp.pm_with_19;

			// reduce the sets of intensities
			new_tp.intens.resize(new_tp.frag_idxs.size());
			for (f=0; f<new_tp.frag_idxs.size(); f++)
			{
				new_tp.intens[f].clear();
				if (org_tp.intens[f].size()>0)
				{
					new_tp.intens[f].resize(selected_length+1,NEG_INF);
					int j;
					for (j=0; j<selected_length; j++)
					{
						new_tp.intens[f][j]=org_tp.intens[f][j+best_idx];
					}

				}
			}

			if (verbose)
			{
				cout << "\t";
				new_tp.print(config);
				cout << endl;
			}
		}
		if (verbose)
			cout << endl;
	}
}