spectrum.cpp

MS-Clustering is designed to rapidly cluster large MS/MS datasets. The program merges similar spectr

	}

	if (peaks.size() >0)
	{
		min_peak_mass = peaks[0].mass -1; // margin 1 Dalton
		max_peak_mass = peaks[peaks.size()-1].mass +1;
	}
	else
		return false;


	if (charge>=2)
	{
		max_mass =  org_pm_with_19 > max_peak_mass ? 
					org_pm_with_19 : max_peak_mass;
	}
	else
		max_mass = max_peak_mass;

	max_mass += 11.0; // margin of error

	size_idx = config->calc_size_idx(charge,org_pm_with_19);

	return true;
}

bool Spectrum::init_from_QCPeaks(Config* _config, void *QCPeaks_ptr, int num_peaks, 
								 void *ssf_ptr, bool ind_perfrom_filtering)
{
	SingleSpectrumFile *ssf = (SingleSpectrumFile *)ssf_ptr;
	QCPeak     *QCpeaks = (QCPeak *)QCPeaks_ptr;

	peaks.clear();
	config = _config;
	file_name = "xxx";

	charge = (ssf->charge>0 ? ssf->charge : 2);
	m_over_z = ssf->m_over_z;
	org_pm_with_19 = m_over_z * charge - ( charge - 1)* 1.0025;
	peptide=ssf->peptide;

	peaks.resize(num_peaks);
	int i;
	for (i=0; i<num_peaks; i++)
	{
		peaks[i].mass = QCpeaks[i].mass;
		peaks[i].intensity = QCpeaks[i].intensity;
		total_original_intensity+=peaks[i].intensity;
	}

	if (peaks.size() >0)
	{
		min_peak_mass = peaks[0].mass -1; // margin 1 Dalton
		max_peak_mass = peaks[peaks.size()-1].mass +1;
	}
	else
		return false;


	if (charge>=2)
	{
		max_mass =  org_pm_with_19 > max_peak_mass ? 
					org_pm_with_19 : max_peak_mass;
	}
	else
		max_mass = max_peak_mass;

	max_mass += 11.0; // margin of error

	size_idx = 0;

	init_spectrum(ind_perfrom_filtering);

	return true;
}


void Spectrum::print_spectrum(ostream& os) const
{
	int i;
	if (file_name.length() > 0)
		os << "#FILE " << file_name << endl;
	if (peptide.get_num_aas()>0)
		os << "#SEQ " << peptide.as_string(config) << endl;

	os << fixed << setprecision(2) << org_pm_with_19 << " " << charge << endl;

	for (i=0; i<peaks.size(); i++)
	{
		os << left << setw(5) << i << setw(8) << setprecision(NUM_SIG_DIGITS) 
			<< fixed << right  << peaks[i].mass;
		os << setw(12) << right  << setprecision(1) << peaks[i].intensity << " " << setw(4) 
			<< setprecision(1) << peaks[i].iso_level <<  setw(4) 
			<< setprecision(1) << peaks[i].log_local_rank << "\t" << setprecision(3) << exp(peaks[i].log_rand_prob) << endl;
	}
}


void Spectrum::output_as_MGF(ostream& os) const
{
	os << "BEGIN IONS" << endl;
	os << "TITLE=" << file_name << endl;
	
	if (peptide.get_num_aas()>0)
		os << "SEQ=" << peptide.as_string(config) << endl;
	
	if (scan_number>=0 && cluster_size == 1)
		os << "SCAN=" << scan_number << endl;
	
//	if (retention_time >=0)
//		os << "RT=" << retention_time << endl;

	if (cluster_size>1)
		os << "CLUSTER_SIZE=" << cluster_size << endl;

	os << "CHARGE=+" << charge << endl;

	os << "PEPMASS=" << fixed << setprecision(NUM_SIG_DIGITS) << m_over_z << endl;

	int i;
	for (i=0; i<peaks.size(); i++)
		os << fixed << setprecision(NUM_SIG_DIGITS) << peaks[i].mass << " " 
		   << fixed << setprecision(NUM_SIG_DIGITS) << peaks[i].intensity << endl;
	
	os << "END IONS" << endl << endl;
}



void Spectrum::print_expected_by(ostream& os) const
{
	vector<string> labels;
	labels.push_back("b");
	labels.push_back("y");
	print_expected_fragment_peaks(labels,os);
}

void Spectrum::print_expected_fragment_peaks(vector<string>& frag_labels, ostream& os) const
{
	int i;
	const vector<FragmentType>& all_fragments = config->get_all_fragments();
	vector<mass_t> break_masses;
	
	const int frag_label_width =20;

	if (peptide.get_num_aas() == 0)
		return;

	vector<string> pre_frags, suf_frags;
	vector<int> pre_frag_idxs, suf_frag_idxs;
	
	for (i=0; i<frag_labels.size(); i++)
	{
		if (frag_labels[i].length()>0 &&
			(frag_labels[i][0] == 'a' || frag_labels[i][0]== 'b' || frag_labels[i][0] == 'c') )
		{
			pre_frags.push_back(frag_labels[i]);
			pre_frag_idxs.push_back(config->get_frag_idx_from_label(frag_labels[i]));
		}
		else
		{
			suf_frags.push_back(frag_labels[i]);
			suf_frag_idxs.push_back(config->get_frag_idx_from_label(frag_labels[i]));
		}
	}


	peptide.calc_expected_breakage_masses(config,break_masses);
	mass_t true_mass_with_19 = peptide.get_mass() + MASS_OHHH;
	os << peptide.as_string(config) << " (" << fixed << setprecision(4) << true_mass_with_19 << ")" << endl;

	if ( suf_frags.size() == 0)
	{
		os << setw(frag_label_width*frag_labels.size()+3)<< setfill('-') << right << " " << endl;
		os << setfill(' ');

		os << "|  |";
		for (i=0; i<frag_labels.size(); i++)
		{
			int w = (frag_label_width+6- frag_labels[i].length())/2;
			os << setw(w) << " " << setw(frag_label_width-w-1) << left << frag_labels[i] << "|";
		}
		os<<endl;

		os << setw(frag_label_width*frag_labels.size()+3)<< setfill('-') << right << " " << endl;
		os << setfill(' ');

		
		for (i=0; i<break_masses.size(); i++)
		{
			int region = config->calc_region_idx(break_masses[i],
				true_mass_with_19, charge, min_peak_mass, max_peak_mass);

			const RegionalFragments& rf = config->get_regional_fragments(charge,size_idx,region);

			if (rf.get_frag_type_idxs().size() == 0)
			{
				cout << "Error: no fragments selected for region " << charge << " " <<
						size_idx << " " << region << endl;
				exit(1);
			}

			cout << "|" << setw(2) << left << i <<"|";
			
			int j;
			for (j=0; j<frag_labels.size(); j++)
			{
				int idx=pre_frag_idxs[j];
				if (rf.get_position_of_frag_type_idx(idx)<0)
					idx=-1;
				
				if (idx<0)
				{
					os << setw(frag_label_width-2) <<  left << "      - ";
				}
				else
				{
					mass_t exp_mass = all_fragments[idx].calc_expected_mass(break_masses[i],
										true_mass_with_19);
					int p_idx=get_max_inten_peak(exp_mass,config->get_tolerance());
					os << " " << setw(6) << right << fixed << setw(9) << setprecision(3) << exp_mass;
					if (p_idx>=0)
					{
						os << " " << setw(6) << fixed << setprecision(3) << peaks[p_idx].mass - exp_mass;
					}
					else
						os << "       ";
				}
				os << " |";
			}
			os << endl;
		}
		os << setw(frag_label_width*frag_labels.size()+5)<< setfill('-') << right << " " << endl;
		os << setfill(' ');
	}
	else
	if ( pre_frags.size() == 0)
	{
		os << setw(frag_label_width*frag_labels.size()+5)<< setfill('-') << right << " " << endl;
		os << setfill(' ');

		os << "|  |";
		for (i=0; i<frag_labels.size(); i++)
		{
			int w = (frag_label_width - 1- frag_labels[i].length())/2;
			os << setw(w) << " " << setw(frag_label_width-1-w) << left << frag_labels[i] << "|";
		}
		os<<endl;

		os << setw(frag_label_width*frag_labels.size()+5)<< setfill('-') << right << " " << endl;
		os << setfill(' ');

		for (i=0; i<break_masses.size(); i++)
		{
			int region = config->calc_region_idx(break_masses[i],
				true_mass_with_19, charge, min_peak_mass, max_peak_mass);

			const RegionalFragments& rf = config->get_regional_fragments(charge,size_idx,region);

			cout << "|" << setw(2) << left << break_masses.size() - i - 1<<"|";
			
			int j;
			for (j=0; j<frag_labels.size(); j++)
			{
				int idx=suf_frag_idxs[j];
				if (rf.get_position_of_frag_type_idx(idx)<0)
					idx=-1;

				if (idx<0)
				{
					os << setw(13) <<  left << "      - ";
				}
				else
				{
					mass_t exp_mass = all_fragments[idx].calc_expected_mass(break_masses[i],
										true_mass_with_19);
					int p_idx=get_max_inten_peak(exp_mass,config->get_tolerance());
					os << " " << setw(6) << right << fixed << setw(9) << setprecision(3) << exp_mass;
					if (p_idx>=0)
					{
						os << " " << setw(6) << fixed << setprecision(3) << peaks[p_idx].mass - exp_mass;
					}
					else
						os << "       ";
				}
				os << " |";
			}
			os << endl;
		}

		os << setw(frag_label_width*frag_labels.size()+5)<< setfill('-') << right << " " << endl;
		os << setfill(' ');
	}
	else // have both prefix and suffix fragments, separate between them
	{
		os << setw(frag_label_width*frag_labels.size()+7)<< setfill('-') << right << " " << endl;
		os << setfill(' ');

		os << "|  |";
		for (i=0; i<pre_frags.size(); i++)
		{
			int w = (frag_label_width - 2 - pre_frags[i].length())/2;
			os << setw(w) << " " << setw(frag_label_width-2-w) << left << pre_frags[i] << "|";
		}
		os << "   |";
		for (i=0; i<suf_frags.size(); i++)
		{
			int w = (frag_label_width - 2 - suf_frags[i].length())/2;
			os << setw(w) << " " << setw(frag_label_width-2-w) << left << suf_frags[i] << "|";
		}
		os<<endl;

		os << setw(frag_label_width*frag_labels.size()+7)<< setfill('-') << right << " " << endl;
		os << setfill(' ');

		for (i=0; i<break_masses.size(); i++)
		{
			int region = config->calc_region_idx(break_masses[i],
				true_mass_with_19, charge, min_peak_mass, max_peak_mass);

			const RegionalFragments& rf = config->get_regional_fragments(charge,size_idx,region);


			cout << "|" << setw(2) << left << i <<"|";
			
			int j;
			for (j=0; j<pre_frags.size(); j++)
			{
				int idx=pre_frag_idxs[j];
				if (rf.get_position_of_frag_type_idx(idx)<0)
					idx=-1;

				if (idx<0)
				{
					os << setw(13) <<  left << "      - ";
				}
				else
				{
					mass_t exp_mass = all_fragments[idx].calc_expected_mass(break_masses[i],
										true_mass_with_19);
					int p_idx=get_max_inten_peak(exp_mass,config->get_tolerance());
					os << " " << setw(6) << right << fixed << setw(9) << setprecision(3) << exp_mass;
					if (p_idx>=0)
					{
						os << " " << setw(6) << fixed << setprecision(3) << peaks[p_idx].mass - exp_mass;
					}
					else
						os << "       ";
				}
				os << " |";
			}

			// output suffix fragments

			cout << " " << setw(2) << left << break_masses.size() - i - 1<<"|";
			for (j=0; j<suf_frags.size(); j++)
			{
				int idx=suf_frag_idxs[j];
				if (rf.get_position_of_frag_type_idx(idx)<0)
					idx=-1;

				if (idx<0)
				{
					os << setw(13) <<  left << "      - ";
				}
				else
				{
					mass_t exp_mass = all_fragments[idx].calc_expected_mass(break_masses[i],
										true_mass_with_19);
					int p_idx=get_max_inten_peak(exp_mass,config->get_tolerance());
					os << " " << setw(6) << right << fixed << setw(9) << setprecision(3) << exp_mass;
					if (p_idx>=0)
					{
						os << " " << setw(6) << fixed << setprecision(3) << peaks[p_idx].mass - exp_mass;
					}
					else
						os << "       ";
				}
				os << " |";
			}


			os << endl;
		}

		os << setw(frag_label_width*frag_labels.size()+7)<< setfill('-') << right << " " << endl;
		os << setfill(' ');
	}

	// print pm stats

	

	os << "true - org_pm:    " << fixed << setprecision(4) << true_mass_with_19 << " - " <<
		org_pm_with_19 << "  = " << true_mass_with_19 - org_pm_with_19 << endl;
	
	if (this->corrected_pm_with_19>0)
	{
		os << "true - cor1_pm:   "  << fixed << setprecision(4) << true_mass_with_19 << " - " <<
		corrected_pm_with_19 << "  = " << true_mass_with_19 - corrected_pm_with_19 << endl;
	}
	if (this->secondary_pm_with_19>0)
	{
		os << "true - cor2_pm: "  << fixed << setprecision(4) << true_mass_with_19 << " - " <<
		secondary_pm_with_19 << "  = " << true_mass_with_19 - secondary_pm_with_19 << endl;
	}
	os << endl;
	
}


// checks several charges to see which one has a good m over z match
// writes the expected b/y ions if finds a match
bool Spectrum::check_m_over_z_and_sequence(ostream& os)
{
	int c;
	mass_t pep_mass = peptide.get_mass() + MASS_OHHH;

	for (c=1; c<=4; c++)
	{
		mass_t pm_19 = m_over_z * c - (c-1)*MASS_PROTON;

	//	cout <<"c: " <<c << "   pm19: " << pm_19 << endl;

		if (fabs(pm_19-pep_mass)<7)
		{
			set_org_pm_with_19(pm_19);
			set_charge(c);
		//	print_expected_by(os);
		//	os << endl;
			return true;
		}
	}

	for (c=1; c<=4; c++)
	{
		mass_t pm_19 = m_over_z * c - (c-1)*MASS_PROTON;

		cout <<"c: " <<c << "   pm19: " << pm_19 << endl;
	}

	os << "No charge matches: " << peptide.as_string(config) << "  (" <<
		peptide.get_mass() + 19.083 << ") !!!! " << endl << endl;

	return false;
}


ostream& operator << (ostream& os, const Spectrum& spec)
{
	const vector<Peak>& peaks = spec.get_peaks();
	const string& file_name = spec.get_file_name();
	const Peptide& peptide = spec.get_peptide();
	const Config *config = spec.get_config();

	if (file_name.length() > 0)
		os << "#FILE " << file_name << endl;
	if (peptide.get_num_aas()>0)
		os << "#SEQ " << peptide.as_string(config) << endl;

	os << spec.get_org_pm_with_19() << " " << spec.get_charge() << endl;
	int i;
	for (i=0; i<peaks.size(); i++)
		os << setw(8) << left << peaks[i].mass << " \t" << peaks[i].intensity << 
		" \t" << exp(peaks[i].log_rand_prob)<< endl;

	return os;
}