fragmentation.cpp

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

#include "Fragmentation.h"
#include "Config.h"


const char *neutral_loss_labels[]={"","-H2O","-NH3","-H2OH2O","-H2ONH3","-NH3NH3"};
const mass_t neutral_loss_offsets[]={0,-MASS_H2O,-MASS_NH3,-MASS_H2OH2O,-MASS_H2ONH3,-MASS_NH3NH3};
const int num_neutral_losses = sizeof(neutral_loss_labels)/sizeof(char *);

const char  *prefix_base_labels[]={"a","b","c"};
const mass_t prefix_base_offsets[]={-26.9871,MASS_PROTON,18.0343}; 
const int num_prefix_bases = sizeof(prefix_base_labels)/sizeof(char *);

const char  *suffix_base_labels[]={"x","y","z"};
const mass_t suffix_base_offsets[]={45.9976,MASS_OHHH,2.99};
const int num_suffix_bases = sizeof(suffix_base_labels)/sizeof(char *);


void FragmentType::read_fragment(istream& is)
{
	char line[128];
	is.getline(line,128);
	istringstream iss(line);
	char direction;
	iss >> direction;
	if (direction == 'p')
	{
		orientation = PREFIX;
	}
	else if (direction == 's')
	{
		orientation = SUFFIX;
	}
	else
	{
		cout << "Error reading fragment: " << line << endl;
		exit(1);
	}
	iss >> charge >> offset >> label;
}

void FragmentType::write_fragment(ostream& os) const
{
	if (orientation == PREFIX)
	{
		os << "p ";
	}
	else
		os << "s ";

	os << charge << " ";
	os << fixed << setprecision(5) << offset << " ";
	os << label << endl;
}

/********************************************************
Creates a string label from the fragment information.
*********************************************************/
void FragmentType::make_frag_label(mass_t tolerance)
{
	

	label = (orientation == PREFIX ? "p" : "s");

	if (charge>1)
	{
		ostringstream os;
		os << charge;
		label+= os.str();
	}

	if (offset != 0)
	{
		ostringstream os;
		os << fixed << setprecision(1) << offset;
		if (offset>0)
			label += "+";
		label += os.str();
	}

	// use standard labels for known fragments


	char **base_labels;
	mass_t *base_offsets;
	int num_bases;

	if (orientation == PREFIX)
	{
		base_labels  = (char **)prefix_base_labels;
		base_offsets = (mass_t *)prefix_base_offsets;
		num_bases    = num_prefix_bases;
	}
	else
	{
		base_labels  = (char **)suffix_base_labels;
		base_offsets = (mass_t *)suffix_base_offsets;
		num_bases    = num_suffix_bases;
	}

	// check all possible fragments
	mass_t tol = tolerance;
	
	if (tolerance>=0.5)
	{
		tol = (0.8 * tolerance)/charge;
	}
	else if (tolerance>0.1)
		tol = tolerance / charge;
	
	int b;
	for (b=0; b<num_bases; b++)
	{
		int n;
		for (n=0; n<num_neutral_losses; n++)
		{
			mass_t calc_offset = (base_offsets[b] + (charge - 1)*MASS_PROTON + neutral_loss_offsets[n])/(mass_t)charge;
			if (fabs(calc_offset-offset)<tol)
			{
				label = base_labels[b];
				if (charge>1)
					label+=char('0'+charge);
				
				if (n>0)
					label+=neutral_loss_labels[n];

				offset = (base_offsets[b] + neutral_loss_offsets[n] + (charge - 1)*MASS_PROTON) / charge;
				break;
			}
		}
		if (n<num_neutral_losses)
			break;
	}
}	
	



// removes fragments that appear to be isotopic peaks of previously
// selected fragments such as b+1, y+2 etc.
void FragmentTypeSet::remove_isotopic_fragments(mass_t tolerance)
{

	sort_fragments_according_to_probs();

	int f;
	for (f=0; f<this->fragments.size(); f++)
	{
		FragmentType& frag = fragments[f];

		if (frag.label.length()<1)
			frag.make_frag_label(tolerance);

		// don't check known labels
		if (frag.label[0] != 's' &&  frag.label[0] != 'p')
			continue;

		int j;
		for (j=0; j<f; j++)
		{
			FragmentType& previous = fragments[j];

			if (previous.orientation == frag.orientation &&
				previous.charge      == frag.charge)
			{
				mass_t offset_diff = frag.offset - previous.offset;

				offset_diff *= frag.charge;

				if (offset_diff<3.3) // we do not want unrecognized fragments with such
								     // a small distance from selected fragments
				{
					cout << " --- removing isotopic frag " << frag.label << endl;
					frag.prob = -1;
					frag.spec_count = 0;
					break;
				}
			}
		}
	}

	sort_fragments_according_to_probs();

	while (fragments.size()>0 && fragments[fragments.size()-1].prob<0)
		fragments.pop_back();

}


void FragmentTypeSet::sort_fragments_according_to_probs()
{
	sort(fragments.begin(),fragments.end());
}

/*************************************************
Outputs labels in single line
**************************************************/
void FragmentTypeSet::output_fragment_labels(ostream& os) const
{
	int i;
	for (i=0; i<fragments.size()-1; i++)
		os << fragments[i].label << " ";
	os << fragments[i].label << endl;
}


void FragmentTypeSet::print() const
{
	int i;
	
	cout << setw(4) << left << " "
		     << setw(10) << left << "Label " 
		//	 << setw(5) << "Pidx"
		//	 <<  "offset"
			 << endl;
	for (i=0; i<fragments.size(); i++)
		cout << setw(4) << left <<i 
		     << setw(10) << left << fragments[i].label
		//	 << setw(5) << fragments[i].parent_frag_idx
		//	 << fragments[i].offset_from_parent_frag
			 << endl;
}






/***************************************************************
For each fragment idx its parent fragment is defined as the
fragment with the highest probability that has the same charge
and oreientation as the fragment in question.
****************************************************************/
void FragmentTypeSet::set_parent_frag_idxs()
{
	int i;

	for (i=0; i<fragments.size(); i++)
	{
		int j;
		for (j=0; j<i; j++)
		{
			if (fragments[j].charge == fragments[i].charge &&
				fragments[j].orientation == fragments[i].orientation)
				break;
		}

		// found parent_frag_idx
		if (j<i)
		{
			fragments[i].parent_frag_idx = j;
			fragments[i].offset_from_parent_frag = fragments[i].offset - fragments[j].offset;
		}
		else
		{
			fragments[i].parent_frag_idx=-1;
			fragments[i].offset_from_parent_frag =0;
		}
	}
}






void RegionalFragments::select_fragments_with_minimum_prob(score_t min_prob, int max_num_frags)
{
	while (frag_type_idxs.size()>0 && frag_probs[frag_type_idxs.size()-1]<min_prob)
	{
		frag_type_idxs.pop_back();
		frag_probs.pop_back();
	}

	if (max_num_frags>0)
	{
		while (frag_type_idxs.size()>max_num_frags)
		{
			frag_type_idxs.pop_back();
			frag_probs.pop_back();
		}
	}

}


void FragmentCombo::print_combo(const Config *config, ostream & os) const
{
	int i;
	for (i=0; i<this->frag_inten_idxs.size(); i++)
		os << config->get_fragment(frag_inten_idxs[i]).label << " ";
	os << "|";
	for (i=0; i<this->frag_no_inten_idxs.size(); i++)
		os << " " << config->get_fragment(frag_no_inten_idxs[i]).label;
	os << endl;
}


void FragmentCombo::read_combo(Config *config, istream& is)
{
	frag_inten_idxs.clear();
	frag_no_inten_idxs.clear();

	vector<string> tokens;
	char buff[128];
	is.getline(buff,128);

	istringstream iss(buff);
	string s;
	while (iss >> s)
		tokens.push_back(s);

	// find pos of separator
	int i;
	for (i=0; i<tokens.size(); i++)
		if (! strcmp(tokens[i].c_str(),"|"))
			break;
	if (i== tokens.size())
	{
		cout << "Error: couldn't find '|' for FragmentCombo : " << buff << endl;
		exit(1);
	}
	const int sep_pos = i;
	
	// get inten idxs
	for (i=0; i<sep_pos; i++)
	{
		int f_idx=config->get_frag_idx_from_label(tokens[i]);
		if (f_idx>=0)
		{
			frag_inten_idxs.push_back(f_idx);
		}
		else
			break;
	}
	
	// get no inten idxs
	for (i=sep_pos+1; i<tokens.size(); i++)
	{
		int f_idx=config->get_frag_idx_from_label(tokens[i]);
		if (f_idx>=0)
		{