quickclusteringspectra.cpp

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

#include "QuickClustering.h"
#include "auxfun.h"
#include "AnnotatedSpectrum.h"


// static member dclr

vector<QCPeak> ClusterSpectrum::tmp_peak_area1;
vector<QCPeak> ClusterSpectrum::tmp_peak_area2;

vector<int> ClusterSpectrum::min_num_occurences;
int   ClusterSpectrum::num_top_peaks_per_1000_da;
float ClusterSpectrum::large_num_occurence_ratios;
int   ClusterSpectrum::large_cluster_size;


void ClusterSpectrum::filter_peaks_with_slidinig_window()
{
	vector<bool> inds;
	vector<QCPeak> new_peaks;

	mark_top_peaks_with_sliding_window(&peaks[0],
									   peaks.size(), 
									   config->get_local_window_size(),
									   config->get_max_number_peaks_per_local_window(),
									   inds);

	int i;
	for (i=0; i<inds.size(); i++)
	{
		if (inds[i])
			new_peaks.push_back(peaks[i]);
	}

	peaks=new_peaks;

}


struct peak_idx_pair
{
	bool operator< (const peak_idx_pair& other) const
	{
		return (intensity>other.intensity);
	}
	int idx;
	intensity_t intensity;
};
/***********************************************************************
Uses a heuristic approach jumps every half window
************************************************************************/
bool mark_top_peaks_with_sliding_window(const QCPeak *peaks, 
										int num_peaks, 
										mass_t window_size, 
										int num_peaks_per_window, 
										vector<bool>& indicators)
{
	// filter low intensity noise
	// and mark those that are good peaks
	const mass_t half_window_size = 0.5 * window_size;
	const int max_peak_idx = num_peaks -1;

	if (num_peaks<=5)
	{
		indicators.resize(num_peaks,true);
		return false;
	}
	int i;
	for (i=0; i<5; i++)
	{
		if (peaks[i].scaled_intensity<=0)
			break;
	}

	const bool use_scaled_intensity = (i==5);
	int start_window_idx =0;

	indicators.resize(num_peaks,false);
	indicators[0]=true;
	indicators[max_peak_idx]=true;

	while (start_window_idx<max_peak_idx)
	{
		const mass_t max_window_mass = peaks[start_window_idx].mass + window_size;

		int end_window_idx=start_window_idx;
		while (end_window_idx<max_peak_idx && peaks[end_window_idx].mass<max_window_mass)
			end_window_idx++;


		if (end_window_idx - start_window_idx>num_peaks_per_window)
		{
			const int num_peaks_in_window = end_window_idx - start_window_idx+1;
			vector<peak_idx_pair> pairs;
			pairs.resize(num_peaks_in_window);

			if (use_scaled_intensity)
			{
				int i;
				for (i=0; i<num_peaks_in_window; i++)
				{
					const int peak_idx = i+start_window_idx;
					peak_idx_pair& pair = pairs[i];
					pair.idx = peak_idx ;
					pair.intensity = peaks[peak_idx].scaled_intensity;
				}
			}
			else
			{
				int i;
				for (i=0; i<num_peaks_in_window; i++)
				{
					const int peak_idx = i+start_window_idx;
					peak_idx_pair& pair = pairs[i];
					pair.idx = peak_idx ;
					pair.intensity = peaks[peak_idx].intensity;
				}
			}

			sort(pairs.begin(),pairs.end());

			if (pairs[0].intensity<pairs[1].intensity)
			{
				printf("Error: with peak intensity order (possible corruption in the files)!\n");
			//	int i;
			//	for (i=0; i<pairs.size(); i++)
			//		cout << i << " " << pairs[i].intensity << endl;
			//	exit(1);
				return false;
			}

			int i;
			for (i=0; i<num_peaks_per_window; i++)
				indicators[pairs[i].idx]=true;	
		}
		else 
		{
			int i;
			for (i=start_window_idx; i<=end_window_idx; i++)
				indicators[i]=true;
		}

		// advance half a window
		const mass_t mid_mass = peaks[start_window_idx].mass + half_window_size;
		start_window_idx++;
		while (start_window_idx<max_peak_idx && peaks[start_window_idx].mass<mid_mass)
			start_window_idx++;

	}

	return true;
}






/****************************************************************
*****************************************************************/
void ClusterSpectrum::create_new_cluster(Config *config,
										 BasicSpectrum& bs,
										 int cluster_idx)
{
	int i;
	this->tmp_cluster_idx = cluster_idx;
	this->config = config;
	this->tolerance = config->get_tolerance();
	this->m_over_z = bs.ssf->m_over_z;
	this->peptide_str = bs.ssf->peptide.as_string(config);

	this->num_spectra_in_cluster =1;

	if (bs.ssf->sqs >=0)
		this->best_sqs_spec_idx = 0;

	retention_time = bs.ssf->retention_time;

	maximum_good_peaks_to_output = (int)( (config->get_number_of_strong_peaks_per_local_window()/
									 config->get_local_window_size()) 
									 * bs.peaks[bs.num_peaks-1].mass);

	maximum_peaks_vector_size = (int)(2.5*maximum_good_peaks_to_output);

	if (bs.num_peaks>maximum_peaks_vector_size)
	{
		maximum_peaks_vector_size = bs.num_peaks;
	}

	peaks.reserve(maximum_peaks_vector_size);
	peaks.resize(bs.num_peaks);
	for (i=0; i<bs.num_peaks; i++)
		peaks[i]=bs.peaks[i];

	bs.ssf->assigned_cluster = cluster_idx;
	basic_spectra.push_back(bs);	
}






/*************************************************************************
// adds the spectrum to this cluster.
**************************************************************************/
void ClusterSpectrum::add_spectrum_to_cluster(BasicSpectrum& bs, 
											  const vector<int>& spec_top_idxs,
											  float top_x_masses[NUM_TOP_CLUSTER_PEAKS])
{
	int i;

	bs.ssf->assigned_cluster = tmp_cluster_idx;
	basic_spectra.push_back(bs);
	const mass_t tolerance = (config->get_tolerance()>0.1? config->get_tolerance()*0.8 : config->get_tolerance());

	if (bs.ssf->sqs<0 || basic_spectra.size()>MAX_SIZE_FOR_SQS_REP)
	{
		add_peak_list(bs.peaks,bs.num_peaks,tolerance,1);
		static vector<int> cluster_top_idxs;
		set_adjusted_inten(&peaks[0],peaks.size());
		set_cluster_m_over_z();
		select_top_peak_idxs(&peaks[0],peaks.size(),m_over_z,tolerance,
				cluster_top_idxs, top_peak_masses, num_top_peaks_per_1000_da, config);
		set_top_ranked_idxs(cluster_top_idxs);
	}
	else
	{
		if (basic_spectra.size() == MAX_SIZE_FOR_SQS_REP)
		{
			create_consensus_by_binning_basic_spectra();
		}
		else
		{
			if (best_sqs_spec_idx<0)
			{
				cout << "Error: best_sqs_idx<0 ! " << endl;
				best_sqs_spec_idx=0;
			}

			if (bs.ssf->sqs>basic_spectra[best_sqs_spec_idx].ssf->sqs)
			{
				best_sqs_spec_idx = basic_spectra.size()-1;
				peaks.resize(bs.num_peaks);
				int i;
				for (i=0; i<bs.num_peaks; i++)
					peaks[i]=bs.peaks[i];
				
				m_over_z = bs.ssf->m_over_z;
				set_top_ranked_idxs(spec_top_idxs);
				set_top_masses(top_x_masses);
			}
		}
	}



	/*	// update retention time
	if (retention_time>0)
	{
		retention_time =0;
		int count=0;
		for (i=0; i<basic_spectra.size(); i++)
		{
			const float &rt = basic_spectra[i].ssf->retention_time;
			if (rt>=0)
			{
				count++;
				retention_time += rt;
			}
		}
		retention_time /= (float)count;
	}*/

	// look for a consensus string
	vector<string> peps;
	vector<int> counts;

	for (i=0; i<basic_spectra.size(); i++)
	{
		if (basic_spectra[i].ssf->peptide.get_num_aas()>3)
		{
			string pep_str = basic_spectra[i].ssf->peptide.as_string(config);
			int j;
			for (j=0; j<peps.size(); j++)
				if (! strcmp(pep_str.c_str(),peps[j].c_str()) )
					break;
				
			if (j==peps.size())
			{
				peps.push_back(pep_str);
				counts.push_back(1);
			}
			else
				counts[j]++;
		}
	}

	if (peps.size() == 1)
	{
		this->peptide_str = peps[0];
	}
}


/*************************************************************************
 tries to add the cluster
 succeeds only if the similarity of the two originals to the new consensus
 is above the sim_tresh (returns true if it made the addition, false otherwise)
**************************************************************************/
bool ClusterSpectrum::add_cluster(ClusterSpectrum& cs, float sim_thresh)
{
	const int size_before = basic_spectra.size();

	int i;
	for (i=0; i<cs.basic_spectra.size(); i++)
	{
		BasicSpectrum& bs = cs.basic_spectra[i];

		bs.ssf->assigned_cluster = tmp_cluster_idx;
		basic_spectra.push_back(bs);
	}

	mass_t tolerance = (config->get_tolerance()>0.1? config->get_tolerance()*0.8 : config->get_tolerance());


	if (basic_spectra.size()<=MAX_SIZE_FOR_SQS_REP)
	{
		if (best_sqs_spec_idx<0)
		{
		//	cout << "Error: found best_sqs_idx<0 !!" << endl;
			best_sqs_spec_idx = 0;
		}

		float best_sqs = basic_spectra[best_sqs_spec_idx].ssf->sqs;
		int new_idx=-1;
		int i;
		for (i=size_before; i<basic_spectra.size(); i++)
		{
			if (basic_spectra[i].ssf->sqs > best_sqs)
			{
				new_idx = i;
				best_sqs = basic_spectra[i].ssf->sqs;
			}
		}

		if (new_idx>0)
		{
			best_sqs_spec_idx = new_idx;
			peaks.resize(cs.peaks.size());
			int i;
			for (i=0; i<cs.peaks.size(); i++)
				peaks[i]=cs.peaks[i];

			m_over_z = cs.get_m_over_z();
			set_top_masses(cs.get_top_peak_masses());
			set_top_ranked_idxs(cs.get_top_ranked_idxs());
		}
	}
	else
	{
		// avoid creating consensus for very large clusters
		if (size_before <= MAX_SIZE_FOR_SQS_REP && 
			cs.basic_spectra.size() <= MAX_SIZE_FOR_SQS_REP)
		{
			create_consensus_by_binning_basic_spectra(); 
		}
		else
		{
			add_peak_list(cs.get_peaks_pointer(),
						  cs.get_num_peaks(),tolerance, 
					      cs.num_spectra_in_cluster);

			static vector<int> cluster_top_idxs;
			set_adjusted_inten(&peaks[0],peaks.size());

			set_cluster_m_over_z();

			select_top_peak_idxs(&peaks[0],peaks.size(),m_over_z,tolerance,
				cluster_top_idxs, top_peak_masses, num_top_peaks_per_1000_da, config);

			set_top_ranked_idxs(cluster_top_idxs);
		}
	}



	// look for a consensus string
	vector<string> peps;
	vector<int> counts;

	for (i=0; i<basic_spectra.size(); i++)
	{
		if (basic_spectra[i].ssf->peptide.get_num_aas()>3)
		{
			string pep_str = basic_spectra[i].ssf->peptide.as_string(config);
			int j;
			for (j=0; j<peps.size(); j++)
				if (! strcmp(pep_str.c_str(),peps[j].c_str()) )
					break;
			
			if (j==peps.size())
			{
				peps.push_back(pep_str);
				counts.push_back(1);
			}
			else
				counts[j]++;

		}
	}

	if (peps.size() == 1)
	{
		this->peptide_str = peps[0];
	}

	return true;
}

/*************************************************************************
// adds the given list to the clusters current list.
// if the new cluster is larger than X, peak scores are weighted according
// to the percentage in which the peaks appear
// list of peaks is then filtered to remove excess peaks
**************************************************************************/
bool ClusterSpectrum::add_peak_list(
		const QCPeak *second_peaks, 
		int num_second_peaks, 
		mass_t tolerance, 
		int num_basic_spectra_added,
		bool need_to_scale)
{
	static vector<QCPeak> tmp_peaks1, tmp_peaks2;
	static vector<float> scaling_factors;

	if (scaling_factors.size()<100)
	{
		scaling_factors.resize(101);
		int i;
		for (i=0; i<=100; i++)
			scaling_factors[i]=0.2 + 0.2 *pow(1+i*0.01,5);
	}


	num_spectra_in_cluster+=num_basic_spectra_added;

	const int    num_joined = this->peaks.size() + num_second_peaks;	

	if (tmp_peaks1.size() < num_joined)
	{
		tmp_peaks1.resize(2*num_joined);
		tmp_peaks2.resize(2*num_joined);
	}

	if (num_joined<=0)
		return true;

	// place a merged list of peaks in tmp_peaks1
	const int num_peaks = peaks.size();
	int a_idx=0, b_idx=0, p_idx=0;

	while (a_idx<num_peaks && b_idx<num_second_peaks)
	{
		if (peaks[a_idx].mass<= second_peaks[b_idx].mass)
		{
			tmp_peaks1[p_idx++]=peaks[a_idx++];
		}
		else
			tmp_peaks1[p_idx++]=second_peaks[b_idx++];
	}

	while (a_idx<num_peaks)
		tmp_peaks1[p_idx++]=peaks[a_idx++];

	while (b_idx<num_second_peaks)
		tmp_peaks1[p_idx++]=second_peaks[b_idx++];

	// use 3 rounds with increasing tolerance to join peaks
	// each time put them in the new list in area2
	vector<mass_t> tolerances;
	tolerances.resize(3,0);
	tolerances[0]=tolerance*0.25;
	tolerances[1]=tolerance*0.5;
	tolerances[2]=tolerance;
	
	int round;
	int num_area1 = p_idx;
	for (round=0; round<3; round++)
	{
		const mass_t join_tolerance = tolerances[round];
		vector<QCPeak>& area1 = ( (round == 1) ? tmp_peaks2 : tmp_peaks1);
		vector<QCPeak>& area2 = ( (round == 1) ? tmp_peaks1 : tmp_peaks2);
		int j=0;
		int i;

		area2[0]=area1[0];
		for (i=1; i<num_area1; i++)
		{
			if (area1[i].mass - area2[j].mass<join_tolerance)
			{
				const QCPeak& peak1 = area1[i];
				QCPeak& peak2 = area2[j];

				const int num_occurences = peak2.num_occurences + peak1.num_occurences;

				intensity_t sum_intens = peak1.intensity + peak2.intensity;
				mass_t weight = peak1.intensity / sum_intens;
				peak2.mass = (peak1.mass * weight) + (1.0-weight)* peak2.mass;
				peak2.intensity = sum_intens;
			
				peak2.num_occurences = num_occurences;
			}
			else
				area2[++j]=area1[i];
		}
		num_area1 = j+1;
	}

	const int merged_num_peaks = num_area1;

	// scale intensity according to peak probability
	if (need_to_scale && num_spectra_in_cluster>2)
	{
		const float idx_mult = 100.0 / num_spectra_in_cluster;
		int i;

		for (i=0; i<merged_num_peaks; i++)
		{
			QCPeak& peak = tmp_peaks2[i];
			int idx = (int)(idx_mult*peak.num_occurences);
			if (idx>100)
				idx=100;

			peak.scaled_intensity = peak.intensity * scaling_factors[idx];
		}
	}

	// check if we can just copy the peaks
	if (merged_num_peaks <= maximum_peaks_vector_size)
	{
		// copy over peaks
		peaks.clear();
		int i;
		for (i=0; i<merged_num_peaks; i++)
			peaks.push_back(tmp_peaks2[i]);

		return true;
	}

	vector<bool> indicators;
	mark_top_peaks_with_sliding_window(
		&tmp_peaks2[0],
		merged_num_peaks,
		config->get_local_window_size(),
		(int)(config->get_max_number_peaks_per_local_window()*2.5),
		indicators);

//	cout << config->get_local_window_size() << "\t" <<  config->get_max_number_peaks_per_local_window() << endl;

	peaks.clear();
	int i;
	for (i=0; i<merged_num_peaks; i++)
		if (indicators[i])
			peaks.push_back(tmp_peaks2[i]);

	return true;
}


/***********************************************************************
// recursively merges the peak lists from the various spectra
// performs a merge until the pointers list has only one entry
************************************************************************/
void ClusterSpectrum::merge_peak_lists(vector<QCPeak>& org_peaks,
									   vector<QCPeak>& new_peaks,
									   vector<PeakListPointer>& pointers)
{
	QCPeak *org_peak_area = &org_peaks[0];
	QCPeak *new_peak_area = &new_peaks[0];
	while (pointers.size()>1)
	{
		vector<PeakListPointer> new_pointers;
		QCPeak *n_pos = new_peak_area;
		new_pointers.clear();