quickclusteringspectra.cpp

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

			sim_sums[i]+=sim;
			sim_sums[j]+=sim;
		}
	}

	float max_sim=-1;
	int max_sim_idx=-1;
	for (i=0; i<basic_spectra.size(); i++)
	{
		if (sim_sums[i]>max_sim)
		{
			max_sim=sim_sums[i];
			max_sim_idx=i;
		}
	}

	
	peaks.resize(basic_spectra[max_sim_idx].num_peaks);
	for (i=0; i<basic_spectra[max_sim_idx].num_peaks; i++)
		peaks[i]=basic_spectra[max_sim_idx].peaks[i];

	set_top_ranked_idxs(top_idxs[max_sim_idx]);
	set_cluster_m_over_z();

	return max_sim_idx;
}



int ClusterSpectrum::select_max_sqs_spectrum_as_consensus()
{
	int max_idx=0;
	float max_sqs=0;
	int i;

	for (i=0; i<basic_spectra.size(); i++)
	{
		if (basic_spectra[i].ssf->sqs>max_sqs)
		{
			max_idx=i;
			max_sqs = basic_spectra[i].ssf->sqs;
		}
	}

	BasicSpectrum& spec = basic_spectra[max_idx];
	float top_x_masses[NUM_TOP_CLUSTER_PEAKS];
	vector<int> top_idxs;

	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, top_x_masses, 20);

	peaks.resize(basic_spectra[max_idx].num_peaks);
	for (i=0; i<basic_spectra[max_idx].num_peaks; i++)
		peaks[i]=basic_spectra[max_idx].peaks[i];


	set_top_ranked_idxs(top_idxs);
	set_cluster_m_over_z();

	return max_idx;
}

/************************************************************************
// sets the m_over_z as the average of the m_over_z of the basic_spectra
*************************************************************************/
void ClusterSpectrum::set_cluster_m_over_z()
{
	int i;

	m_over_z=0;

/*	int c = basic_spectra[0].ssf->charge;
	for (i=1; i<basic_spectra.size(); i++)
		if (basic_spectra[i].ssf->charge != c)
			break;

	// not all charges are the same, use the first m_over_z 
	if (i<basic_spectra.size())
	{
		m_over_z = basic_spectra[0].ssf->m_over_z;
		return;
	} */

	for (i=0; i<basic_spectra.size(); i++)
		m_over_z += basic_spectra[i].ssf->m_over_z;

	if (basic_spectra.size() == 0)
	{
		cout << "Error: cluster has no basic spectra!" << endl;
		exit(1);
	}
	m_over_z /= basic_spectra.size();
}


/************************************************************************
Choose Majority charge
*************************************************************************/
void ClusterSpectrum::set_charge()
{

	vector<int> charge_counts;
	charge_counts.resize(8,0);

	if (basic_spectra.size() == 1)
	{
		charge = basic_spectra[0].ssf->charge;
		return;
	}

	int i;
	for (i=0; i<basic_spectra.size(); i++)
		charge_counts[basic_spectra[i].ssf->charge]++;
	
	int max_charge=0;
	for (i=0; i<8; i++)
		if (charge_counts[i]>=charge_counts[max_charge])
			max_charge=i;

	charge = max_charge;

	if (charge ==0)
		return;

	const float ratio_max = (float)charge_counts[max_charge]/(float)basic_spectra.size();

	if (charge>1 && ratio_max<=0.5)
		charge=0;

	if (charge == 2 || charge == 3)
	{
		float ratio2 = (float)charge_counts[2]/(float)basic_spectra.size();
		float ratio3 = (float)charge_counts[3]/(float)basic_spectra.size();

		if (ratio2>=0.33 && ratio3>=0.33)
			charge=0;
	}
}



// creates the title (file name) for the cluster
void ClusterSpectrum::make_title(string& name, int batch_idx, int cluster_idx)
{
	ostringstream os1;
	ostringstream os2;

	os1 << batch_idx;
	os2 << cluster_idx;

	title = name + "." + os1.str() + "." + os2.str();
}


// finds how many spectra have a peptide sequence that 
// doesn't match the majority assignment
// adds spectra to the mismatch category only if they contain 
// masses that are not within tolerance of each other (modulo -3,-2,-1,+1,+2,+3)
int  ClusterSpectrum::get_num_misassigned_spectra(mass_t pm_tolerance, int* mismatched_with_pep) const
{
	vector<string> peptide_strings;
	vector<int> counts;
	int i;

	if (mismatched_with_pep)
		*mismatched_with_pep=0;

	for (i=0; i<basic_spectra.size(); i++)
	{
		if (basic_spectra[i].ssf->peptide.get_num_aas()<2)
			continue;

		const string& pep_str = basic_spectra[i].ssf->peptide.as_string(config);
	
		int j;
		for (j=0; j<peptide_strings.size(); j++)
		{
			if (! strcmp(peptide_strings[j].c_str(),pep_str.c_str()))
			{
				counts[j]++;
				break;
			}
		}

		if (j<peptide_strings.size())
			continue;

		peptide_strings.push_back(pep_str);
		counts.push_back(1);
	}

	if (counts.size()<=1)
		return 0;

	int max_idx=0;
	int max_count=counts[0];

	for (i=1; i<peptide_strings.size(); i++)
	{
		if (counts[i]>max_count)
		{
			max_idx=i;
			max_count=counts[i];
		}
	}

	Peptide max_pep;
	string max_pep_str = peptide_strings[max_idx];
	max_pep.parse_from_string(config,max_pep_str);
	mass_t max_pep_mass_with_19 = max_pep.get_mass() + MASS_OHHH;

	int num_missmatched=0;
	int num_missmatched_with_peptide=0;
	for (i=0; i<basic_spectra.size(); i++)
	{
		string pep_str = basic_spectra[i].ssf->peptide.as_string(config);
		if (pep_str.length()>0 && ! strcmp(pep_str.c_str(),max_pep_str.c_str()) )
			continue;

		// check that if parent mass is wrong
		mass_t pep_m_over_z = basic_spectra[i].ssf->m_over_z;
		mass_t best_off = 9999;
		mass_t best_pep_mass = -1;
		int c;
		for (c=1; c<=4; c++) // assume the spectrum's charge is between 1-4
		{
			mass_t pep_mass = pep_m_over_z * c - (c-1)* MASS_PROTON;
			mass_t pm_offset;
			for (pm_offset=-4.0; pm_offset<=4.0; pm_offset+=1.0)
			{
				mass_t offset = fabs(max_pep_mass_with_19 - pep_mass - pm_offset);
				if (offset<best_off)
				{
					best_off = offset;
					best_pep_mass = pep_mass;
				}
			}
		}

		if (best_off>pm_tolerance)
		{
		//	cout << setprecision(8) << max_pep_mass_with_19 << " " << best_pep_mass << "  " << 
		//		max_pep_mass_with_19 - best_pep_mass << endl;
			num_missmatched++;

			if (pep_str.length()>0)
			{
			//	cout << pep_str << " " << num_missmatched_with_peptide << endl;
				num_missmatched_with_peptide++;
			}
		}
	}

	if (mismatched_with_pep)
		*mismatched_with_pep=num_missmatched_with_peptide;

	return num_missmatched;
}



// returns true if there is a mjority annotation (with 75% of the annotated spectra)
// if so, returns its string and mass in the variables
bool ClusterSpectrum::has_majority_annotation(string& consensus_pep_str, mass_t& consensus_pep_mass) const
{
	vector<string> peptide_strings;
	vector<int> counts;
	int num_annoatated_spectra=0;
	int i;

	peptide_strings.clear();
	counts.clear();
	

	for (i=0; i<basic_spectra.size(); i++)
	{
		if (basic_spectra[i].ssf->peptide.get_num_aas()<2)
			continue;

		const string& pep_str = basic_spectra[i].ssf->peptide.as_string(config);
		num_annoatated_spectra++;
	
		int j;
		for (j=0; j<peptide_strings.size(); j++)
		{
			if (! strcmp(peptide_strings[j].c_str(),pep_str.c_str()))
			{
				counts[j]++;
				break;
			}
		}

		if (j<peptide_strings.size())
			continue;

		peptide_strings.push_back(pep_str);
		counts.push_back(1);
	}

	if (counts.size()==0)
		return false;

//	cout << ">>> " << basic_spectra.size() << " <" << counts.size() << "> " <<
//		num_annoatated_spectra << endl;

	

	int max_idx=0;
	int max_count=counts[0];

	for (i=1; i<counts.size(); i++)
	{
		if (counts[i]>max_count)
		{
			max_idx=i;
			max_count=counts[i];
		}
	}

//	cout <<  max_count << " " << num_annoatated_spectra << endl;
	if ( ((double)max_count/(double)num_annoatated_spectra) <0.75)
		return false;

	Peptide max_pep;

	consensus_pep_str = peptide_strings[max_idx];
	max_pep.parse_from_string(config,consensus_pep_str);
	consensus_pep_mass = max_pep.get_mass() + MASS_OHHH;
	
	return true;
}


int ClusterSpectrum::get_num_basic_spectra_with_peptide() const
{
	int num_peps=0;
	int i;

	for (i=0; i<basic_spectra.size(); i++)
		if (basic_spectra[i].ssf->peptide.get_num_aas()>1)
			num_peps++;
	
	return num_peps;
}

	
void ClusterSpectrum::print_cluster_peptides() const
{
	vector<string> peptide_strings;
	vector<int> counts;
	int i;

	for (i=0; i<basic_spectra.size(); i++)
	{
		if (basic_spectra[i].ssf->peptide.get_num_aas()<3)
			continue;

		const string& pep_str = basic_spectra[i].ssf->peptide.as_string(config);

		int j;
		for (j=0; j<peptide_strings.size(); j++)
		{
			if (! strcmp(peptide_strings[j].c_str(),pep_str.c_str()))
			{
				counts[j]++;
				break;
			}
		}

		if (j<peptide_strings.size())
			continue;

		peptide_strings.push_back(pep_str);
		counts.push_back(1);
	}

	
	int num_with_peptide=0;
	for (i=0; i<counts.size(); i++)
		num_with_peptide+=counts[i];

	cout << "Cluster: " << this->tmp_cluster_idx << "  " << basic_spectra.size() << " (";
	cout << num_with_peptide << ")" << endl;
	for (i=0; i<counts.size(); i++)
		cout << setw(4) << left << counts[i] << peptide_strings[i].c_str() << endl;
}



void ClusterSpectrum::print_cluster_similarities()
{
	int i;
	mass_t tolerance = config->get_tolerance();

	for (i=0; i<basic_spectra.size(); i++)
	{
		BasicSpectrum& spec = basic_spectra[i];
		float top_x_masses[NUM_TOP_CLUSTER_PEAKS];
		vector<int> spec_top_idxs;
		set_adjusted_inten(spec.peaks,spec.num_peaks);
		select_top_peak_idxs(spec.peaks,spec.num_peaks,spec.ssf->m_over_z,
			tolerance,spec_top_idxs, top_x_masses, num_top_peaks_per_1000_da);


		float sim = calc_selected_dot_prod(tolerance,
					&peaks[0], peaks.size(), top_ranked_peak_idxs,
					spec.peaks,spec.num_peaks, spec_top_idxs);

		cout << i << " " << sim << " " << spec.ssf->single_name;
	}
}



void BasicSpectrum::output_to_mgf(ostream& mgf, Config *config, const char *seq) const
{
	mgf << "BEGIN IONS" << endl;
	mgf << "TITLE=" <<  ssf->single_name << endl;
	
	if (ssf->peptide.get_num_aas()>0)
	{
		mgf << "SEQ=" << ssf->peptide.as_string(config) << endl;
	}
	else if (seq && strlen(seq)>2)
		mgf << "SEQ=" << seq << endl;
	
	if (ssf->type == MZXML)
	{
//		MZXML_single *mzxml_single = (MZXML_single *)ssf;
//		if (mzxml_single->scan_number>=0)
//			mgf << "SCAN=" <<mzxml_single->scan_number << endl;

//		if (mzxml_single->retention_time>=0)
//			mgf << "RT=" << mzxml_single->retention_time << endl;
	}

	mgf << "CHARGE=+" << ssf->charge << endl;
		
	mgf << "PEPMASS=" << ssf->m_over_z << endl;
	
	int i;
	for (i=0; i<this->num_peaks; i++)
		mgf << fixed << setprecision(3) << peaks[i].mass << " " << peaks[i].intensity << endl;

	mgf << "END IONS" << endl << endl;
}



void BasicSpectrum::output_to_mgf(FILE* mgf_stream, Config *config, const char *seq) const
{
	fprintf(mgf_stream,"BEGIN IONS\n");
	fprintf(mgf_stream,"TITLE=%s\n",ssf->single_name.c_str());
	
	if (ssf->peptide.get_num_aas()>0)
	{
		fprintf(mgf_stream,"SEQ=%s\n",ssf->peptide.as_string(config).c_str());
	}
	else if (seq && strlen(seq)>2)
		fprintf(mgf_stream,"SEQ=%s\n", seq);
	
	if (ssf->type == MZXML)
	{
		MZXML_single *mzxml_single = (MZXML_single *)ssf;
		if (mzxml_single->scan_number>=0)
			fprintf(mgf_stream,"SCAN=%d\n",mzxml_single->scan_number);
	}

	fprintf(mgf_stream,"CHARGE=+%d\n",ssf->charge);
		
	fprintf(mgf_stream,"PEPMASS=%.3f\n",ssf->m_over_z);
	
	int i;
	for (i=0; i<this->num_peaks; i++)
		fprintf(mgf_stream,"%.3f %.3f\n",peaks[i].mass, peaks[i].intensity);

	fprintf(mgf_stream,"END IONS\n\n");
}




void ClusterSpectrum::write_cluster_basic_spectra_to_mgf(char *mgf_file) const
{
	ofstream mgf(mgf_file);
	int i;
	for (i=0; i<basic_spectra.size(); i++)
		basic_spectra[i].output_to_mgf(mgf,config);
	mgf.close();
}

// writes the spectrum to the output file in the mgf format
void ClusterSpectrum::write_spectrum_to_mgf(ostream& mgf,bool write_peak_count) const
{
	mgf << "BEGIN IONS" << endl;
	mgf << "TITLE=" <<  title << endl;

	if (peptide_str.length()>0)
		mgf << "SEQ=" << peptide_str << endl;


//	if (retention_time>=0)
//		mgf << "RT=" << retention_time << endl;
//	mgf << "CLUSTER_SIZE=" << basic_spectra.size() << endl;	

	if (this->charge>0 && config->get_use_spectrum_charge())
		mgf << "CHARGE=" << this->charge << "+" << endl;
		
	mgf << "PEPMASS=" << this->m_over_z << endl;

	// filter peaks if there are too many

	if (this->basic_spectra.size()>1)
	{
		if (peaks.size() <= maximum_good_peaks_to_output)
		{	
			int i;
			for (i=0; i<peaks.size(); i++)
			{
				mgf << fixed << setprecision(3) << peaks[i].mass << " " << setprecision(2) << 
					(peaks[i].scaled_intensity>0 ? peaks[i].scaled_intensity : peaks[i].intensity);
				if (write_peak_count)
					mgf << "   " << peaks[i].num_occurences << "/" << peaks[i].max_num_occurences;
				mgf << endl;
			}
		}
		else
		{
			vector<bool> good_peak_indicators;
			mark_top_peaks_with_sliding_window(&peaks[0], peaks.size(),
				config->get_local_window_size(),
				config->get_max_number_peaks_per_local_window(),
				good_peak_indicators);

			int i;
			for (i=0; i<peaks.size(); i++)
			{
				if (! good_peak_indicators[i])
					continue;

				mgf << fixed << setprecision(3) << peaks[i].mass << " " << setprecision(2) << 
					(peaks[i].scaled_intensity>0 ? peaks[i].scaled_intensity : peaks[i].intensity);

				if (write_peak_count)
					mgf << "   " << peaks[i].num_occurences << "/" << peaks[i].max_num_occurences;
				mgf << endl;
			}
		}
	}
	else // output all peaks in the single basic spectrum - no filtering
	{
		const QCPeak *peaks = basic_spectra[0].peaks;
		const int num_peaks = basic_spectra[0].num_peaks;
		int i;
		for (i=0; i<num_peaks; i++)
				mgf << fixed << setprecision(3) << peaks[i].mass << " " << setprecision(3) << 
					peaks[i].intensity << endl;
	}

	mgf << "END IONS" << endl << endl;
	
}


// writes the spectrum to the output file in the mgf format
void ClusterSpectrum::write_spectrum_to_pkl_single(string& file_name) const
{
	ofstream pkl(file_name.c_str());
	if (! pkl.is_open())
	{
		cout << "Error: couldn't open pkl for writing " << file_name << endl;
		exit(1);
	}

	double total_ion_current =0;
	int i;
	for (i=0; i<basic_spectra.size(); i++)
		total_ion_current += basic_spectra[i].ssf->precursor_intensity;

//	cout << "TIC: " << total_ion_current << endl;

	pkl << setprecision(4) << fixed << m_over_z << "\t";
	pkl << setprecision(0) << total_ion_current << "\t";