annotatedspecturm.cpp

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

		if (peak_annotations[i].size()>0)
			n++;

	return n;
}


float AnnotatedSpectrum::get_explianed_intensity() const
{
	if (peak_annotations.size()==0 || (peaks.size() != peak_annotations.size()) )
	{
		cout << "Error: mismatch in peaks annotations!" << endl;
		exit(1);
	}

	float tot_inten=0;
	float ann_inten=0;
	int i;
	for (i=0; i<peaks.size(); i++)
	{
		tot_inten += peaks[i].intensity;
		if (peak_annotations[i].size()>0)
			ann_inten += peaks[i].intensity;
	}
	
	return (ann_inten/tot_inten);
}

/**************************************************************************
Checks if the spectrum has a suffcient stretch of b or y ladders
***************************************************************************/
bool AnnotatedSpectrum::has_stretch_of_b_or_y(int min_stretch_length, int max_skip)
{
	vector<mass_t> b_ladder,y_ladder;
	int b_frag_idx = config->get_frag_idx_from_label("b");
	int y_frag_idx = config->get_frag_idx_from_label("y");
	const vector<mass_t>& aa2mass = config->get_aa2mass();
	mass_t tolerance = config->get_tolerance();

	if (breakages.size() ==0)
	{
		cout << "Error: need to first annotate spectrum!" << endl;
		exit(1);
	}

	b_ladder.resize(breakages.size(),0);
	b_ladder[0]=1;
	b_ladder[b_ladder.size()-1]=get_true_mass()+1;
	y_ladder.resize(breakages.size(),0);
	y_ladder[0]=get_true_mass_with_19();
	y_ladder[y_ladder.size()-1]=19;

	int i;
	vector<mass_t> exp_b_masses,exp_y_masses;
	peptide.calc_expected_breakage_masses(config,exp_b_masses);
	exp_y_masses = exp_b_masses;

	for (i=0; i<exp_y_masses.size(); i++)
		exp_y_masses[i] = get_true_mass_with_19() - exp_y_masses[i];
	for (i=0; i<exp_b_masses.size(); i++)
		exp_b_masses[i]+=0.95;


	for (i=0; i<breakages.size(); i++)
	{
		int b_pos= breakages[i].get_position_of_frag_idx(b_frag_idx);
		if (b_pos>=0)
			b_ladder[i]=breakages[i].fragments[b_pos].mass;

		int y_pos = breakages[i].get_position_of_frag_idx(y_frag_idx);
		if (y_pos>=0)
			y_ladder[i]=breakages[i].fragments[y_pos].mass;
	}

	
	const vector<int>& amino_acids = peptide.get_amino_acids();

//	for (i=0; i<b_ladder.size(); i++)
//		cout << y_ladder[i] << " " << exp_y_masses[i] << endl;


	// check for stretch
	for (i=0; i<= b_ladder.size()-min_stretch_length; i++)
	{
		int gaps=0;
		int j;

		mass_t last_b = exp_b_masses[i];
		mass_t exp_mass =0;
		for (j=0; j<min_stretch_length; j++)
		{
			if (j>0)
				exp_mass += aa2mass[amino_acids[i+j-1]];

			if (b_ladder[i+j]>0 && (j==0 || 
									fabs(b_ladder[i+j]-last_b - exp_mass)<=0.5 * tolerance ) )
			{
				gaps=0;
				exp_mass =0;
				last_b = exp_b_masses[i+j];
			}
			else
				gaps++;

			if (gaps>max_skip)
				break;
		}

		if (j==min_stretch_length)
			return true;
	} 

	// check for stretch
	int num_aa = amino_acids.size();
	for (i=0; i<= y_ladder.size()-min_stretch_length; i++)
	{
		int gaps=0;
		int j;
		mass_t last_y = exp_y_masses[i];
		mass_t exp_mass =0;
		for (j=0; j<min_stretch_length; j++)
		{
			if (j>0)
				exp_mass += aa2mass[amino_acids[i+j-1]];

			if (y_ladder[i+j]>0 && (j==0 || 
									fabs(last_y - y_ladder[i+j]- exp_mass)<=0.5*tolerance ) )
			{
				gaps=0;
				exp_mass =0;
				last_y = exp_y_masses[i+j];
			}
			else
				gaps++;

			if (gaps>max_skip)
				break;
		}
		if (j==min_stretch_length)
			return true;
	}

//	exit(0);
	return false;
}


/********************************************************************************
Extracts tables of peak intensities and peak masses.
#rows = all fragments
# column = # num breakages (peptide length + 1)
*********************************************************************************/
void AnnotatedSpectrum::extract_annotated_intens_and_masses(
											 vector< vector<intensity_t> >& intens,
											 vector< vector<mass_t> >&      masses) const
{
	const int num_frags = config->get_all_fragments().size();
	const int num_breakages = breakages.size();
	int i;
	intens.resize(num_frags);
	masses.resize(num_frags);
	for (i=0; i<num_frags; i++)
	{
		intens[i].resize(num_breakages,0);
		masses[i].resize(num_breakages,0);
	}

	for (i=0; i<peak_annotations.size(); i++)
	{
		int j;
		for (j=0; j<peak_annotations[i].size(); j++)
		{
			const PeakAnnotation& ann = peak_annotations[i][j];
			if (ann.breakage_idx>0 && ann.frag_type_idx<num_frags)
			{
				intens[ann.frag_type_idx][ann.breakage_idx]=peaks[i].intensity;
				masses[ann.frag_type_idx][ann.breakage_idx]=peaks[i].mass;
			}
		}
	}
}


// report file names of files that have a minimal b/y stretch
void filter_file_list(const FileManager& fm, Config *config, ostream& os)
{
	FileSet fs;
	fs.select_all_files(fm);
	int good=0;
	while(1)
	{
		AnnotatedSpectrum as;

		if (! fs.get_next_spectrum(fm,config,&as))
			break;

		as.annotate_spectrum(as.get_true_mass_with_19());
	//	as.print_expected_by();
		if (as.has_stretch_of_b_or_y(8,1))
		{
			good++;
			cout << good << " ";
			os << as.get_file_name() << endl;
			
		}
		else
			cout << "XXXX"<<endl;
	}

//	cout << good << " / " << fs.get_total_spectra() << " = " << 
//			good / (double)fs.get_total_spectra() << endl;
}




void AnnotatedSpectrum::print_annotated(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;

	cout << "Peak  Mass        Inten   Rank   Annotation" << endl;
	int strong_idx=0;
	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(2) << peaks[i].intensity << " " << setw(4) 
		//	<< setprecision(1) << peaks[i].iso_level <<  setw(4) 
		//	<< setprecision(1) << peaks[i].log_local_rank 
			<< setw(4) << right << peaks[i].rank;

		if (strong_idx<strong_peak_idxs.size())
		{
			if (strong_peak_idxs[strong_idx]==i)
			{
				strong_idx++;
				os << " *";
			}
			else
				os << "  ";
		}
		else
			os << "  ";

		int j;
		for (j=0; j<peak_annotations[i].size(); j++)
			os << "  " << peak_annotations[i][j].label;
			
		os << endl;
	}
}


void print_dataset_spectra_by_stats(Config *config, char *mgf_file)
{
	FileManager fm;
	FileSet     fs;

	fm.init_from_mgf(config,mgf_file);
	fs.select_all_files(fm);

	int b_frag_idx = config->get_frag_idx_from_label("b");
	int y_frag_idx = config->get_frag_idx_from_label("y");

	int counter=0;
	while (1)
	{
		AnnotatedSpectrum as;

		if (! fs.get_next_spectrum(fm,config,&as))
			break;

		as.annotate_spectrum(as.get_true_mass_with_19());

		if (counter==0)
			as.print_expected_by();

		int num_b = as.get_num_observed_frags(b_frag_idx);
		int num_y = as.get_num_observed_frags(y_frag_idx);
		float exp_int = as.get_explianed_intensity();

		
		cout << counter++ << " " << num_b << " " << num_y << " " << exp_int << endl;
	}
}