quickclusteringspectra.cpp

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

		int i;
		for (i=0; i<pointers.size(); i+=2)
		{
			// merge the two lists
			if (i<pointers.size()-1)
			{
				QCPeak* first_list = pointers[i].peaks;
				QCPeak* second_list = pointers[i+1].peaks;
				QCPeak* merge_start = n_pos;

				const int n1 = pointers[i].num_peaks;
				const int n2 = pointers[i+1].num_peaks;
			
				// merge
				int i1=0,i2=0;
				while (i1<n1 && i2<n2)
				{
					if (first_list[i1].mass < second_list[i2].mass)
					{
						*n_pos++=first_list[i1++];
					}
					else
						*n_pos++=second_list[i2++];
				}

				while (i1<n1)
					*n_pos++=first_list[i1++];

				while (i2<n2)
					*n_pos++=second_list[i2++];

				// add pointer for new merged list
				PeakListPointer new_pointer;
				new_pointer.peaks = merge_start;
				new_pointer.num_peaks = n1+n2;
				new_pointers.push_back(new_pointer);
			}
			else // write peaks directly to new area
			{	
				QCPeak* first_list = pointers[i].peaks;
				QCPeak* start_copy = n_pos;
				const int n1 = pointers[i].num_peaks;
				int i1=0;
				while (i1<n1)
					*n_pos++=first_list[i1++];
		
				PeakListPointer new_pointer;
				new_pointer.peaks = start_copy;
				new_pointer.num_peaks = n1;
				new_pointers.push_back(new_pointer);
				n_pos+= n1;
			}
		}

		pointers = new_pointers;

		// switch between the pointers of the peak storage areas
		QCPeak *tmp = org_peak_area;
		org_peak_area = new_peak_area;
		new_peak_area = tmp;
	}
}



/************************************************************************
// joins adjacent peaks, marks invaldiated peaks by assigning their mass to -1
// then condences list to contain only good peaks.
// also calculates for each peak mass what is the maximal number of peaks
// that could be detected at that range (based on all the spectra's min/max
// peak values).
*************************************************************************/
void ClusterSpectrum::join_merged_peak_lists(
										PeakListPointer& plp,
										PeakListPointer& alt_plp,
										int num_merged_spectra,
										mass_t tolerance)
{
	vector<mass_t> join_tolerances;
	int i,t;
	int mid_idx = plp.num_peaks/2;

	// set maximal number of detected peaks
	vector<bool> used_ind;
	used_ind.resize(num_merged_spectra,false);

	// fill from left to middle
	int p_idx,max_peaks_at_idx=0;
	for (p_idx=0; p_idx<=mid_idx; p_idx++)
	{
		QCPeak& peak = plp.peaks[p_idx];
		const int& spec_idx = peak.source_spec_idx;
		if (! used_ind[spec_idx])
		{
			used_ind[spec_idx]=true;
			max_peaks_at_idx++;
		}
		peak.max_num_occurences = max_peaks_at_idx;
	}

	// fill from right to middle
	for (i=0; i<num_merged_spectra; i++)
		used_ind[i]=false;

	max_peaks_at_idx=0;
	for (p_idx = plp.num_peaks-1; p_idx>mid_idx; p_idx--)
	{
		QCPeak& peak = plp.peaks[p_idx];
		const int& spec_idx = peak.source_spec_idx;
		if (! used_ind[spec_idx])
		{
			used_ind[spec_idx]=true;
			max_peaks_at_idx++;
		}
		peak.max_num_occurences = max_peaks_at_idx;
	}

	// join peaks

	join_tolerances.push_back(tolerance*0.15);
	join_tolerances.push_back(tolerance*0.3);
	join_tolerances.push_back(tolerance*0.5);

	for (t=0; t<join_tolerances.size();t++)
	{
		const mass_t join_tolerance = join_tolerances[t];
		const int max_p_idx = plp.num_peaks;

		int p_idx=0, np_idx=0;

		while (p_idx<max_p_idx)
		{
			QCPeak& curr_peak = plp.peaks[p_idx];
			if (curr_peak.mass<0)
				continue;

			int next_idx = p_idx+1;
			while (next_idx<max_p_idx && plp.peaks[next_idx].mass<0)
				next_idx++;
			
			// try joining the peak with the next one ahead
			while (next_idx < max_p_idx &&
				   plp.peaks[next_idx].mass - curr_peak.mass < join_tolerance)
			{
				QCPeak& next_peak = plp.peaks[next_idx];

				// don't try and join peaks that should stay apart
				if (t>=2 &&
					(curr_peak.num_occurences >= curr_peak.max_num_occurences ||
					 next_peak.num_occurences >= next_peak.max_num_occurences ) )
					break;
					 
				const intensity_t sum_inten = curr_peak.intensity + next_peak.intensity;

				mass_t new_mass = (curr_peak.mass * curr_peak.intensity + next_peak.mass * next_peak.intensity) / 
								  sum_inten;

				curr_peak.mass      = new_mass;
				curr_peak.intensity = sum_inten;
				curr_peak.num_occurences += next_peak.num_occurences;

				if (curr_peak.max_num_occurences<next_peak.max_num_occurences)
					curr_peak.max_num_occurences = next_peak.max_num_occurences;
				
				next_peak.mass = -1;
				next_peak.intensity = -100000000;
				next_peak.num_occurences=0;
				
				while (next_idx<max_p_idx && plp.peaks[next_idx].mass<0)
					next_idx++;
			}

			// copy peak to new area
			alt_plp.peaks[np_idx++] = plp.peaks[p_idx];

			// advance to next peak
			p_idx = next_idx;
		}

		alt_plp.num_peaks = np_idx;

		// switch the plps
		PeakListPointer tmp;
		tmp=plp;
		plp=alt_plp;
		alt_plp=plp;
	}
}



/************************************************************************
// selects the consensus peaks - those that appear more than the expected cutoff
// also takes some of the stronger peaks that didn't make the cutoff
// writes the selected peaks into the peaks of the cluster spectrum
*************************************************************************/
void ClusterSpectrum::select_consensus_peaks(PeakListPointer& plp, 
											 PeakListPointer& alt_plp,
											 int num_org_spectra)
{
	vector<bool> keep_indicators;
	keep_indicators.clear();

	keep_indicators.resize(plp.num_peaks,false);


	if (num_org_spectra<this->large_cluster_size)
	{
		int i;
		int num_saved_peaks=0;
		for (i=0; i<plp.num_peaks; i++)
		{
			QCPeak& peak = plp.peaks[i];

			if (peak.num_occurences>= min_num_occurences[peak.max_num_occurences])
				keep_indicators[i]=true;
		}
			
	}
	else // for large clusters might need to use the ratio to calc min num_occurrences
	{
		int i;
		int num_saved_peaks=0;
		for (i=0; i<plp.num_peaks; i++)
		{
			QCPeak& peak = plp.peaks[i];

			if (peak.max_num_occurences < large_cluster_size)
			{
				if (peak.num_occurences>= min_num_occurences[peak.max_num_occurences])
					keep_indicators[i]=true;
			}
			else
			{
				int min_num_occurences = (int)(peak.max_num_occurences * this->large_num_occurence_ratios +0.5);
				if (peak.num_occurences>= min_num_occurences)
					keep_indicators[i]=true;
			}
		}
	}


	// copy peaks
	int i;
	int alt_idx=0;
	QCPeak* peak_list = alt_plp.peaks;
	
	for (i=0; i<plp.num_peaks; i++)
		if (keep_indicators[i])
			peak_list[alt_idx++]=plp.peaks[i];
			
	// filter low intensity peaks

	const mass_t half_window_size = 0.5 * config->get_local_window_size();
	const int num_peaks_in_window = config->get_max_number_peaks_per_local_window();
	int max_peak_idx = alt_idx -1;
	int min_idx=1;
	int max_idx=1;
	
	peaks.clear();
	peaks.reserve(max_peak_idx);

	peaks.push_back(peak_list[0]);

	// check the rest of the peaks
	for (i=1; i<max_peak_idx; i++)
	{
		const mass_t& peak_mass=peak_list[i].mass;
		mass_t min_mass = peak_list[min_idx].mass;
		mass_t max_mass = peak_list[max_idx].mass;

	
		// advance min/max pointers
		while (peak_mass-min_mass > half_window_size)
			min_mass=peak_list[++min_idx].mass;

		while (max_idx < max_peak_idx && max_mass - peak_mass <= half_window_size)
			max_mass=peak_list[++max_idx].mass;

		if (max_mass - peak_mass > half_window_size)
			max_idx--;

		// if there are less than the maximum number of peaks in the window, keep it.
		if (max_idx-min_idx < num_peaks_in_window)
		{
			peaks.push_back(peak_list[i]);
			continue;
		}

		// check if this is one of the top peaks in the window
		int higher_count=0;
		for (int j=min_idx; j<=max_idx; j++)
			if (peak_list[j].intensity > peak_list[i].intensity)
				higher_count++;

		if (higher_count < num_peaks_in_window)
			peaks.push_back(peak_list[i]);
	}
	peaks.push_back(peak_list[max_peak_idx]);
}




void ClusterSpectrum::create_consensus_sepctrum_from_peak_list_pointers(
					  vector<PeakListPointer>& plp, int total_num_peaks)
{
	int i;

	// the merged list pointer is in plp[0]
	merge_peak_lists(tmp_peak_area1,tmp_peak_area2,plp);

	QCPeak *peak_list = plp[0].peaks;
	if (plp[0].num_peaks != total_num_peaks)
	{
		cout << "Error: mismatch in peak numbers of merged lists: " <<
			total_num_peaks << " vs. " << plp[0].num_peaks << endl;
		exit(1);
	}

	
	// assign the appropriate plp pointers
	PeakListPointer merged_peaks_plp = plp[0];
	PeakListPointer alt_plp;
	QCPeak *peak_area1 = &tmp_peak_area1[0];
	QCPeak *peak_area2 = &tmp_peak_area2[0];

	if (merged_peaks_plp.peaks == peak_area1)
	{
		alt_plp.peaks = peak_area2;
	}
	else if (merged_peaks_plp.peaks == peak_area2)
	{
		alt_plp.peaks = peak_area1;
	}
	else
	{
		cout << "Error: mismatch in peak areas!" << endl;
		exit(0);
	}


	join_merged_peak_lists(merged_peaks_plp, alt_plp, 
						   basic_spectra.size(), config->get_tolerance());


	select_consensus_peaks(merged_peaks_plp, alt_plp, basic_spectra.size());


	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);


	// 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];
	}
}




/************************************************************************
	Creates a single consensus spectrum from the basic spectra.
	First merges all the peak lists into a single (sorted) list.
*************************************************************************/
void ClusterSpectrum::create_cluster_by_binning_basic_spectra()
{
	int i;

	create_new_cluster(config,basic_spectra[0],0);
	
	for (i=1; i<basic_spectra.size(); i++)
	{

		add_peak_list(basic_spectra[i].peaks,
			basic_spectra[i].num_peaks,
			config->get_tolerance()*0.8,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);


	// 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];
	}

/*	increase_tmp_storage_size(total_num_peaks);
	
	// write all peaks from the spectra into a single area
	// and create peak list pointer
	int num_peaks_written=0;
	vector<PeakListPointer> plp;
	plp.resize(basic_spectra.size());
	for (i=0; i<basic_spectra.size(); i++)
	{
		QCPeak *list_start  = &tmp_peak_area1[0] + num_peaks_written;
		const int num_peaks_in_list = basic_spectra[i].num_peaks;
		const QCPeak *org_peaks  = basic_spectra[i].peaks;
		plp[i].peaks = list_start;
		plp[i].num_peaks = num_peaks_in_list;

		int j;
		for (j=0; j<num_peaks_in_list; j++)
		{
			QCPeak& tcp = list_start[j];
			tcp.mass = org_peaks[j].mass;
			tcp.intensity = org_peaks[j].intensity;
			tcp.num_occurences=1;
			tcp.source_spec_idx=i;
		}
		num_peaks_written+=num_peaks_in_list;
	}

	create_consensus_sepctrum_from_peak_list_pointers(plp,total_num_peaks);*/

}


/************************************************************************
Clears the current consensus peak list and creates a consensus spectrum
from all cluster members.
*************************************************************************/
void ClusterSpectrum::create_consensus_by_binning_basic_spectra()
{
	int i;
	for (i=0; i<basic_spectra.size(); i++)
	{
		if (i == best_sqs_spec_idx)
			continue;

		add_peak_list(basic_spectra[i].peaks,basic_spectra[i].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);
	
	best_sqs_spec_idx = -1;
}


/************************************************************************
// sets the consensus spectrum to be the basic spectrum with the maximal 
// similarity to other spectra.
*************************************************************************/
int ClusterSpectrum::select_max_similarity_spectrum_as_consensus()
{	
	mass_t tolerance = config->get_tolerance();
	int i;
	vector<float> sim_sums;
	vector< vector<int> > top_idxs;

	top_idxs.resize(basic_spectra.size());
	sim_sums.resize(basic_spectra.size(),0);

	for (i=0; i<basic_spectra.size(); i++)
	{

		BasicSpectrum& spec = basic_spectra[i];
		float top_x_masses[NUM_TOP_CLUSTER_PEAKS];

		

		set_adjusted_inten(spec.peaks,spec.num_peaks);
		select_top_peak_idxs(spec.peaks,spec.num_peaks,spec.ssf->m_over_z,
			tolerance, top_idxs[i], top_x_masses, 20);
	}


	
	for (i=0; i<basic_spectra.size()-1; i++)
	{
		int j;
		for (j=i+1; j<basic_spectra.size(); j++)
		{
			float sim = calc_selected_dot_prod(tolerance,
				basic_spectra[i].peaks,basic_spectra[i].num_peaks, top_idxs[i],
 				basic_spectra[j].peaks,basic_spectra[j].num_peaks, top_idxs[j]);