quickclustering.cpp

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

#include "QuickClustering.h"
#include "RegularRankModel.h"
#include "PMCSQS.h"
#include "auxfun.h"

// the sim matrix stores the similarity distances computed between clusters
int     num_sim_matrix_spectra = 0;
unsigned char  * sim_matrix = NULL;
unsigned char  * max_sim_addr = NULL;



void print_byte(unsigned char byte)
{
	int i;
	unsigned char mask=0X1;

	for (i=0; i<8; i++)
	{
		cout << ( (byte & mask) ? '1' : '0');
		mask <<= 1;
	}
}


void mark_bit_zero(unsigned char *addr, int position)
{
	const unsigned char ANDmasks[]={254,253,251,247,239,223,191,127};
	int bit_offset = (position & 0X7);
	int byte_offset = (position >> 3);
	
	*(addr+byte_offset) &= ANDmasks[bit_offset];
}


void mark_bit_one(unsigned char *addr, int position)
{
	const unsigned char ORmasks[] ={1,2,4,8,16,32,64,128};
	int bit_offset = (position & 0X7);
	int byte_offset = (position >> 3);

	*(addr+byte_offset) |= ORmasks[bit_offset];
}


// returns full int assumes we are at position 31 in the 4 bytes
int get_matrix_32_bits(unsigned char *row_start, int position)
{
	int cell_off = (position >> 5);
	return (*((int *)row_start+cell_off)); 
//	return 1;
}

int get_matrix_val(unsigned char *row_start, int position)
{
	const unsigned char masks[] ={1,2,4,8,16,32,64,128};
	int bit_offset = (position & 0X7);
	int byte_offset = (position >> 3);
	
	return ((*(row_start+byte_offset) & masks[bit_offset]));
//	return 1;
}


// global debugging variable
int wrongly_filtered_spectra;

/***************************************************************************
	This function creates clusters from a list of files containing spectra
	(possibly different file types).
	The cluster spectra are outputted as mgf files in the output dir (x spectra
	per file). In addition, for each cluster file there is a map file that holds
	the indices (position in list, and idx in file) of the original spectra
	that are part of the cluster.
****************************************************************************/
void cluster_full_dataset(Config *config,
							  char *list_file,
							  const string& out_dir,
							  const string& clust_name,
							  int batch_idx,
							  int specs_per_slice,
							  mass_t min_m_over_z,
							  mass_t max_m_over_z,
							  float  min_similarity, 
							  int min_cluster_size,
							  int max_cluster_size,
							  bool verbose,
							  int  max_small_cluster_size,
							  int  k_value,
							  void *pmcsqs,
							  bool ind_pkl_mode,
							  float filter_prob,
							  int  max_mzxml_idx,
							  char *good_anns_file)
{
	const mass_t pm_tolerance = config->get_pm_tolerance();
	const mass_t double_pm_tolerance = pm_tolerance * 2.0;
	const mass_t additional_pm_tolerance = pm_tolerance * 2.5;
	const mass_t tolerance = (config->get_tolerance() <= 0.01) ? config->get_tolerance() :
							  config->get_tolerance() * 0.75;
	FileManager fm;
	FileSet all_spec_fs;
	vector<QCPeak> basic_peaks;      // bulk allocation

	// for cluster histograms
	const int clust_vals[]={1,2,5,10,20,50,100,200,500};
	const int  num_clust_vals = sizeof(clust_vals)/sizeof(int);
	vector<int> clust_counts;

	clust_counts.resize(num_clust_vals+1,0);

	double avg_cluster_size=0;
	int num_clusters=0;
	int total_spectra_in_clusters = 0;
	int spec_idx=0;
	int next_cluster_idx = 0;

	wrongly_filtered_spectra = 0;

	QCOutputter qco;
	qco.init(clust_name,out_dir, batch_idx, min_m_over_z, max_m_over_z,
			 min_similarity, min_cluster_size);

	ostringstream oss;
	oss << batch_idx;
	string batch_str = oss.str();
	string last_good_mass_name = out_dir + "/" + clust_name + "_" + batch_str + "_lastmass.txt";

	ClusterSpectrum::init_statics(config);

	ClusterSpectrum::set_num_top_peaks_per_1000_da(k_value);

	fm.init_from_list_file(config,list_file,min_m_over_z,max_m_over_z);

	all_spec_fs.select_all_files(fm);

	if (max_mzxml_idx>=0)
	{
		const int org_num = all_spec_fs.get_total_spectra();
		cout << "Filtering for max mzxml idx " << max_mzxml_idx << endl;
		all_spec_fs.filter_dat_spectra_by_mzxml_idx(max_mzxml_idx);
		cout << "#spectra went from " << org_num << " to " << all_spec_fs.get_total_spectra() << endl;
	}

	all_spec_fs.sort_according_to_m_over_z();
	const int total_spectra = all_spec_fs.get_total_spectra();
	const vector<SingleSpectrumFile *>& all_ssf = all_spec_fs.get_ssf_pointers();

	if (all_ssf.size()==0)
	{
		cout <<"Warning: no files were selected for clusering in mass range: " <<
			min_m_over_z << " - " << max_m_over_z << endl;
		return;
	}




	// used for debugging purposes
	map<mzXML_annotation,int> ann_map;
	map<mzXML_annotation,int> *ann_map_ptr=NULL;
	if (good_anns_file)
	{
		read_mzXML_annotations_to_map(good_anns_file,ann_map);
		ann_map_ptr = &ann_map;
	}
	
	int total_spectra_read=0;
	int total_mismatched = 0;

	// set the sizes for the static arrays of peaks, sim matrix and determine
	// the maximal slice size for the clustering
	// vlaues are set according to defaults of:
	// 4M spectra
	// 4M peaks
	// Slice of 20000 spectra per iteration
	// Slice width 3 Da.
	

	float multiplier = (float)specs_per_slice / 25000.0;

	int total_spectra_to_cluster = all_ssf.size();
	int num_total_peaks = (int)(4000000 * multiplier); 
	int slice_size      =   specs_per_slice;
	int sim_matrix_size =   ((specs_per_slice+7) / 8);
	int sim_matrix_bytes = sim_matrix_size * sim_matrix_size * 4 + specs_per_slice;



	cout << "Going to cluster " << all_ssf.size() << " spectra in this run." << endl;
	cout << "Need to allocate following memory: " << endl;
	cout << "peaks            " << right << num_total_peaks << endl;
	cout << "spectrum headers " << right << slice_size << endl;
	cout << "sim matrix       " << right << sim_matrix_bytes << " (bytes) " << endl;
	cout << "Using " << k_value << " peaks per 1000 Da for similarity computations." << endl;

	if (pmcsqs)
		cout << "Filtering with quality threshold: " << filter_prob << endl;

	// adjusting the quality filter probability (to improve the default behaviour)
	// singletons should generally have a higher prob for accepting them 
	// this setting is used only for large clustering jobs

	float adjusted_prob_for_min_size_clusters = filter_prob;

	const mass_t mz_range = all_ssf[all_ssf.size()-1]->m_over_z - all_ssf[0]->m_over_z;
	if (mz_range>0)
	{
		float spectra_density =  (all_ssf.size()/mz_range);
		if (spectra_density>1200.0)
		{
			float mult_val = (spectra_density/1600.0)*3.5;
			if (mult_val>3.5)
				mult_val = 3.5;

			adjusted_prob_for_min_size_clusters = (filter_prob * mult_val);
					
			if (adjusted_prob_for_min_size_clusters > 0.35)
				adjusted_prob_for_min_size_clusters = 0.35;

			if (adjusted_prob_for_min_size_clusters<filter_prob)
				adjusted_prob_for_min_size_clusters=filter_prob;
		}
	}
	

	const float filter_prob_for_min_size_clusters = adjusted_prob_for_min_size_clusters;
	const float filter_prob_for_larger_than_min_clusters = filter_prob;	

//	cout << filter_prob_for_min_size_clusters << " " << 
//			filter_prob_for_larger_than_min_clusters << endl;

	
	basic_peaks.resize(num_total_peaks+20000);

	sim_matrix = new unsigned char [sim_matrix_bytes];
	max_sim_addr = sim_matrix + sim_matrix_bytes;

	if (! sim_matrix)
	{
		cout << "Error: couldn't allocate memory for sim matrix!" << endl;
		exit(1);
	}

	double total_sims =0;

	// collect spectra into two ssf vectors
	while (spec_idx<total_spectra)
	{
		const mass_t max_m_over_z = all_ssf[spec_idx]->m_over_z + double_pm_tolerance;
		const mass_t additional_m_over_z = all_ssf[spec_idx]->m_over_z + additional_pm_tolerance;
		int num_used_peaks=0;
		int num_peaks_first_stage = -1;
		int start_idx = spec_idx;
		int end_stage_one_idx = -1;
		int end_additional_idx = -1;

		int start_spec_idx = spec_idx;
		bool add_additional_spectra = false; // flag if spectra from the next mass range
											 // should be added at the later stage
		int i;

		// add spectra wto the 
		int max_spec_idx = spec_idx + slice_size;
		if (max_spec_idx> total_spectra)
			max_spec_idx = total_spectra;

		while (spec_idx < max_spec_idx && 
			   num_used_peaks<num_total_peaks && 
			   all_ssf[spec_idx]->m_over_z<max_m_over_z)
					num_used_peaks+=all_ssf[spec_idx++]->num_peaks;

		end_stage_one_idx = spec_idx-1;
		num_peaks_first_stage = num_used_peaks;

		if (spec_idx < max_spec_idx)
		{
			add_additional_spectra = true;

			while (spec_idx < max_spec_idx && 
			   num_used_peaks<num_total_peaks && 
			   all_ssf[spec_idx]->m_over_z<additional_m_over_z)
					num_used_peaks+=all_ssf[spec_idx++]->num_peaks;

			end_additional_idx = spec_idx-1;
		}

		FileSet cluster_fs, additional_fs;
		cluster_fs.init_from_another_fs(all_spec_fs,start_idx,end_stage_one_idx);

		if (add_additional_spectra)
			additional_fs.init_from_another_fs(all_spec_fs,end_stage_one_idx+1,end_additional_idx);

		vector<ClusterSpectrum> clusters;
		clusters.clear();

		cout << fixed << setprecision(3) << "Clustering: " << all_ssf[start_idx]->m_over_z << " - " << 
			(add_additional_spectra ? all_ssf[end_additional_idx]->m_over_z : 
					all_ssf[end_stage_one_idx]->m_over_z )  << "  (" <<
			spec_idx - start_spec_idx  << "  spectra,  " << num_used_peaks << " peaks)" << endl;

		int num_in_clusters=0;
		total_spectra_read+= cluster_spec_in_file_set(
										config, 
										fm, 
										cluster_fs, 
										tolerance,
										&basic_peaks[0], 
										clusters, 
										min_similarity, 
										max_small_cluster_size,
										max_cluster_size,
										k_value,
										false,
										pmcsqs,
										filter_prob_for_larger_than_min_clusters,
										ann_map_ptr); 


		// join singletons from the next half of the clustering window

		if (add_additional_spectra && additional_fs.get_total_spectra()>0) 
		{
			int num_added =  add_additional_spectra_to_existing_clusters(config,fm,additional_fs,tolerance,
				&basic_peaks[num_peaks_first_stage], clusters, min_similarity,max_cluster_size,
				pmcsqs,filter_prob_for_larger_than_min_clusters, ann_map_ptr, false);

			total_spectra_read += num_added;
		} 

		// update cluster info
		for (i=0; i<clusters.size(); i++)
		{
			ClusterSpectrum& cluster = clusters[i];
			if (cluster.get_tmp_cluster_idx()<0)
				continue;

			if	(cluster.get_num_basic_spectra()<min_cluster_size)
				continue;

			// check if sqs is high enough
			if (pmcsqs && cluster.get_num_basic_spectra() == min_cluster_size &&
				cluster.get_basic_spectrum(0).ssf->sqs < filter_prob_for_min_size_clusters)
				continue;

			total_spectra_in_clusters += cluster.get_num_basic_spectra();

			cluster.set_charge();
			cluster.set_cluster_m_over_z();
	
			if (ind_pkl_mode)
			{
				qco.output_cluster_spectrum_as_single_pkl(cluster);
			}
			else
				qco.output_cluster_spectrum(cluster);

			if (verbose)
			{
				cout << num_clusters<< " " << cluster.get_num_basic_spectra() << endl;
				cluster.print_cluster_similarities();
				cout << endl;
			}

			const int num_spec_in_cluster = cluster.get_num_basic_spectra();
			if (num_spec_in_cluster>=min_cluster_size)
			{
				num_clusters++;
				avg_cluster_size += num_spec_in_cluster;
			}

			// add counts to histogram
			int j;
			for (j=0; j<num_clust_vals; j++)
				if (num_spec_in_cluster<= clust_vals[j])
					break;
			clust_counts[j]++;

			// check for mixed clusters
			int n_mis = clusters[i].get_num_misassigned_spectra();