fragmentation.cpp

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

			frag_no_inten_idxs.push_back(f_idx);
		}
		else
			break;
	}
}




struct idx_prob {
	bool operator< (const idx_prob& other) const
	{
		return (prob>other.prob);
	}
	int idx;
	score_t prob;
};



void RegionalFragments::sort_by_prob()
{
	vector<idx_prob> ip;

	if (frag_type_idxs.size() != frag_probs.size())
	{
		cout << "Error: probs and frag_idxs mismatch!" << endl;
		exit(1);
	}
	
	ip.resize(frag_type_idxs.size());
	int i;
	for (i=0; i<frag_type_idxs.size(); i++)
	{
		ip[i].idx=frag_type_idxs[i];
		ip[i].prob= frag_probs[i];
	}
	sort(ip.begin(),ip.end());
	for (i=0; i<ip.size(); i++)
	{
		frag_type_idxs[i]=ip[i].idx;
		frag_probs[i]=ip[i].prob;
	}

	strong_frag_type_idxs.clear();
	frag_type_combos.clear();
}


/***************************************************
Sets default fragments for PepNovo charge 2
****************************************************/
void RegionalFragments::init_pepnovo_types(int charge, Config *config)
{
	const char* pep_labels[]={"y","b","a","y-H2O","b-H2O","a-H2O","y-NH3","b-NH3",
                           "a-NH3","y-H2OH2O","b-H2OH2O","y-NH3H2O","b-NH3H2O",
						   "y2","b2","y3","b3","y2-H2O","b2-H2O","y2-NH3","b2-NH3",
						   "y3-H2O","b3-H2O","y3-NH3","b3-NH3"};
	frag_type_idxs.resize(13);
	int i;

	for (i=0; i<13; i++)
	{
		string label(pep_labels[i]);
		frag_type_idxs[i]=config->get_frag_idx_from_label(label);
	}
		
	if (charge>1)
	{
		frag_type_idxs.resize(15);
		for (i=13; i<15; i++)
		{
			string label(pep_labels[i]);
			frag_type_idxs[i]=config->get_frag_idx_from_label(label);
		}
	}

	
	if (charge>2)
	{
		frag_type_idxs.resize(21);
		for (i=15; i<21; i++)
		{
			string label(pep_labels[i]);
			frag_type_idxs[i]=config->get_frag_idx_from_label(label);
		}
	}

//	cout << ">> ";
//	for (i=0; i<frag_type_idxs.size(); i++)
//		cout << frag_type_idxs[i] << " ";
//	cout << endl;

	frag_probs.resize(config->get_all_fragments().size(),0);
}




void RegionalFragments::init_with_all_types(int charge, Config *config)
{
	int i;
	frag_type_idxs.clear();
	const vector<FragmentType>& all_fragments = config->get_all_fragments();
	for (i=0; i<all_fragments.size(); i++)
		if (all_fragments[i].charge<=charge)
			frag_type_idxs.push_back(i);

	frag_probs.resize(all_fragments.size(),0);
}






struct idx_pair {
	bool operator< (const idx_pair& other) const
	{
		return (prob>other.prob);
	}
	int f_idx;
	score_t prob;
};

/******************************************************************
// makes the order of the fragments in the frag_idxs and frag_probs 
// vectors be in descending probability order
*******************************************************************/
void RegionalFragments::sort_according_to_frag_probs()
{
	int i;
	vector<idx_pair> pairs;
	if (frag_type_idxs.size() != frag_probs.size())
	{
		cout << "Error: mismtach in frag_idxs and frag_probs sizes!" << 
			frag_type_idxs.size() << " <=> " << frag_probs.size() << endl;
		exit(1);
	}

	pairs.resize(frag_type_idxs.size());
	for (i=0; i<frag_type_idxs.size(); i++)
	{
		pairs[i].f_idx=frag_type_idxs[i];
		pairs[i].prob =frag_probs[i];
	}
	sort(pairs.begin(),pairs.end());
	for (i=0; i<pairs.size(); i++)
	{
		frag_type_idxs[i]=pairs[i].f_idx;
		frag_probs[i]=pairs[i].prob;
	}
}





/*********************************************************************
// chooses all fragments that have high enough probability to be strong
**********************************************************************/
void RegionalFragments::select_strong_fragments(Config *config,
												score_t min_prob, 
												int max_num_strong)
{
	const vector<FragmentType>& all_fragments = config->get_all_fragments();
	int i;

	bool got_prefix = false;
	bool got_suffix = false;

	strong_frag_type_idxs.clear();
	for (i=0; i<frag_type_idxs.size() && i <max_num_strong; i++)
	{
		if (i>=2 && frag_probs[i]<min_prob)
			continue;

		strong_frag_type_idxs.push_back(frag_type_idxs[i]);

		if (all_fragments[frag_type_idxs[i]].orientation == PREFIX)
			got_prefix=true;
		if (all_fragments[frag_type_idxs[i]].orientation == SUFFIX)
			got_suffix=true;
	}

	// try and add an additional frag type...
	if (! got_prefix)
	{
		for (i=0; i<frag_type_idxs.size(); i++)
		{
			if (all_fragments[frag_type_idxs[i]].orientation == PREFIX &&
				frag_probs[i]>min_prob*0.5)
			{
				strong_frag_type_idxs.push_back(frag_type_idxs[i]);
				break;
			}
		}
	}

	if (! got_suffix)
	{
		for (i=0; i<frag_type_idxs.size(); i++)
		{
			if (all_fragments[frag_type_idxs[i]].orientation == SUFFIX &&
				frag_probs[i]>min_prob*0.5)
			{
				strong_frag_type_idxs.push_back(frag_type_idxs[i]);
				break;
			}
		}
	}
}


void Config::set_all_regional_fragment_relationships()
{
	int charge;
	for (charge=0; charge<regional_fragment_sets.size(); charge++)
	{
		int size_idx;
		for (size_idx=0; size_idx<regional_fragment_sets[charge].size(); size_idx++)
		{
			int region_idx;
			for (region_idx=0; region_idx<regional_fragment_sets[charge][size_idx].size(); region_idx++)
				if (regional_fragment_sets[charge][size_idx][region_idx].get_num_fragments()>0)
				{
					regional_fragment_sets[charge][size_idx][region_idx].set_fragment_relationships(this);

				//	if (charge==2 && size_idx==1 && region_idx ==0)
				//		regional_fragment_sets[charge][size_idx][region_idx].print_fragment_relationships(this);
				}
			
		}
	}
}

void Config::print_all_regional_fragment_relationships() const
{
	int charge;
	for (charge=0; charge<regional_fragment_sets.size(); charge++)
	{
		int size_idx;
		for (size_idx=0; size_idx<regional_fragment_sets[charge].size(); size_idx++)
		{
			int region_idx;
			for (region_idx=0; region_idx<regional_fragment_sets[charge][size_idx].size(); region_idx++)
				if (regional_fragment_sets[charge][size_idx][region_idx].get_num_fragments()>0)
				{
					cout << "MODEL " << charge << " " << size_idx << " " << region_idx << endl;
					regional_fragment_sets[charge][size_idx][region_idx].print_fragment_relationships(this);
				}
			
		}
	}
}

/****************************************************************************
This function sets the idxs of parents of the fragments (i.e., they appear in
a higher rank in the table. We distinguish between two cases for the parents:
same charge orientation, and other.
The function also chooses for each strong fragment a 
*****************************************************************************/
void RegionalFragments::set_fragment_relationships(Config *config)
{
	const vector<FragmentType>& all_fragments = config->get_all_fragments();
	const int num_frags = frag_type_idxs.size();

	mirror_frag_idxs.clear();
	parent_idxs.clear();
	parents_with_same_charge_ori_idxs.clear();
	
	mirror_frag_idxs.resize(num_frags);
	parent_idxs.resize(num_frags);
	parents_with_same_charge_ori_idxs.resize(num_frags);
	
	int i;
	for (i=0; i<num_frags; i++)
	{
		const int curr_frag_idx = frag_type_idxs[i];
		const FragmentType& curr_frag = all_fragments[curr_frag_idx];
		int j;
		for (j=0; j<i; j++)
		{
			const int other_frag_idx = frag_type_idxs[j];
			const FragmentType& other_frag = all_fragments[other_frag_idx];

			parent_idxs[i].push_back(other_frag_idx);

			if (curr_frag.charge == other_frag.charge &&
				curr_frag.orientation == other_frag.orientation)
				parents_with_same_charge_ori_idxs[i].push_back(other_frag_idx);
			
			if (curr_frag.orientation != other_frag.orientation)
			{
				const mass_t sum_offset = (curr_frag.offset * curr_frag.charge) + (other_frag.offset * other_frag.charge);
				const int sum_charge = (curr_frag.charge + other_frag.charge);
				if (fabs(sum_offset - MASS_PROTON*(sum_charge-1)-MASS_OHHH)<0.1) 
				{
					mirror_frag_idxs[i].push_back(other_frag_idx);
					mirror_frag_idxs[j].push_back(curr_frag_idx);
				}
			}
		}	
	}
}

void RegionalFragments::print_fragment_relationships(const Config *config) const
{
	const vector<FragmentType>& all_fragments = config->get_all_fragments();
	const int num_frags = frag_type_idxs.size();
	int i;
	for (i=0; i<num_frags; i++)
	{
		const int curr_frag_idx = frag_type_idxs[i];
		const FragmentType& curr_frag = all_fragments[curr_frag_idx];
		int j;

		cout << i <<"\t" << curr_frag.label <<"\tparents  " << parent_idxs[i].size() << " ";
		for (j=0; j<parent_idxs[i].size(); j++)
			cout << "\t" << all_fragments[parent_idxs[i][j]].label;
		cout << endl;
		
		cout <<  "\t\tparents with same " << parents_with_same_charge_ori_idxs[i].size() << " ";
		for (j=0; j<parents_with_same_charge_ori_idxs[i].size(); j++)
			cout << "\t" << all_fragments[parents_with_same_charge_ori_idxs[i][j]].label;
		cout << endl;
		
		cout << "\t\tmirror frags " << mirror_frag_idxs[i].size() << " ";
		for (j=0; j<mirror_frag_idxs[i].size(); j++)
			cout << "\t" << all_fragments[mirror_frag_idxs[i][j]].label;
		cout << endl;

		cout << endl << endl;
	}
}