pmcsqs_model.cpp

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

in the precursor mass.
*****************************************************************************/
int PMCSQS_Scorer::get_optimal_bin(int true_mz_bin, int charge) const
{
	const int max_bin_offset = 6-charge; // look in the range +- of this value
	const vector<PMCRankStats>& pmc_stats = curr_spec_rank_pmc_tables[charge];
	const int min_bin_idx = (true_mz_bin - max_bin_offset>=0 ? true_mz_bin - max_bin_offset : 0);
	const int max_bin_idx = (true_mz_bin + max_bin_offset>= pmc_stats.size() ? pmc_stats.size()-1 :
								true_mz_bin + max_bin_offset);

	if (pmc_stats[true_mz_bin].num_frag_pairs==0 &&
		pmc_stats[true_mz_bin].num_c2_frag_pairs==0)
		return true_mz_bin;
	
	int   optimal_bin_idx=NEG_INF;
	float max_num_pairs=NEG_INF;
	float best_offset=POS_INF;

	if (pmc_stats[true_mz_bin].num_frag_pairs>=pmc_stats[true_mz_bin].num_c2_frag_pairs)
	{
		float max_num_pairs=0;
		int bin_idx;
		for (bin_idx = min_bin_idx; bin_idx<=max_bin_idx; bin_idx++)
			if (pmc_stats[bin_idx].num_frag_pairs > max_num_pairs)
				max_num_pairs = pmc_stats[bin_idx].num_frag_pairs;

		// find minimal offset
		for (bin_idx = min_bin_idx; bin_idx<=max_bin_idx; bin_idx++)
			if (pmc_stats[bin_idx].num_frag_pairs == max_num_pairs &&
				pmc_stats[bin_idx].mean_offset_pairs < best_offset)
			{
				optimal_bin_idx = bin_idx;
				best_offset = pmc_stats[bin_idx].mean_offset_pairs;
			}

		return optimal_bin_idx;
		
	}
	else
	// use the charge 2 fragment pairs
	{
		float max_num_pairs=0; 
		int bin_idx;
		for (bin_idx = min_bin_idx; bin_idx<=max_bin_idx; bin_idx++)
			if (pmc_stats[bin_idx].num_c2_frag_pairs > max_num_pairs)
				max_num_pairs = pmc_stats[bin_idx].num_c2_frag_pairs;

		// find minimal offset
		for (bin_idx = min_bin_idx; bin_idx<=max_bin_idx; bin_idx++)
			if (pmc_stats[bin_idx].num_c2_frag_pairs == max_num_pairs &&
				pmc_stats[bin_idx].mean_offset_c2_pairs < best_offset)
			{
				optimal_bin_idx = bin_idx;
				best_offset = pmc_stats[bin_idx].mean_offset_c2_pairs;
			}

		return optimal_bin_idx;	
	}


	return -1;
}






/*********************************************************************************
Takes a set of samples around the correct mass ([-3+5] every 0.1 Da.)
Selects the bin of the correct mass as positive and a set from offseted m/z
as negative samples. 
**********************************************************************************/
void PMCSQS_Scorer::select_training_sample_idxs(
		int charge,
		const vector<RankBoostSample>& spec_samples,
		const BasicSpectrum& bs,
		int& correct_idx,
		vector<int>& bad_pmc_idxs) const
{
	const vector<PMCRankStats>& pmc_stats = curr_spec_rank_pmc_tables[charge];

	bs.ssf->peptide.calc_mass(config);
	const mass_t pep_mass = bs.ssf->peptide.get_mass()+MASS_H2O;
	const int size_idx = this->get_rank_model_size_idx(charge,pep_mass);
	const mass_t true_mz = (pep_mass + charge)/(float)charge + this->pmc_charge_mz_biases[charge][size_idx];
	const mass_t observed_mz = bs.ssf->m_over_z;

	// check that the training sample has an ok offset
	if (fabs(true_mz-observed_mz)>10.0)
	{
		
		cout << "Erorr in m/z offsets (remove this spectrum from training set): " << endl;
		cout << fixed << setprecision(2) << "file m/z: " << observed_mz << "\t" << 
			"\"true\" m/z: " << true_mz << "\t peptide: " << bs.ssf->peptide.as_string(config) << endl;
		cout << "spectrum: " << bs.ssf->single_name << endl;
		
		cout << "Mass Cys = " << this->config->get_aa2mass()[Cys] << endl;

		exit(1);
	}

	// find the entry with the correct m/z

	int idx=0;
	while (idx<pmc_stats.size() && pmc_stats[idx].m_over_z<true_mz)
		idx++;

	if (idx>= pmc_stats.size())
		idx--;

	if (idx>0 && pmc_stats[idx].m_over_z-true_mz>true_mz-pmc_stats[idx-1].m_over_z)
		idx--;

//	correct_idx=get_optimal_bin(idx,charge);
	correct_idx = idx;

	vector<int> idxs;
	idxs.clear();
	bad_pmc_idxs.clear();


	idxs.push_back(correct_idx+4);
	idxs.push_back(correct_idx+5);
	idxs.push_back(correct_idx+7);
	idxs.push_back(correct_idx+9);
	idxs.push_back(correct_idx+10);
	idxs.push_back(correct_idx+15);

	idxs.push_back(correct_idx+19);
	idxs.push_back(correct_idx+20);
	
	idxs.push_back(correct_idx-4);
	idxs.push_back(correct_idx-5);
	idxs.push_back(correct_idx-7);
	idxs.push_back(correct_idx-9);
	idxs.push_back(correct_idx-10);;
	idxs.push_back(correct_idx-15);
	idxs.push_back(correct_idx-19);
	idxs.push_back(correct_idx-20);

	




	// select upto 5 random samples (make sure they are not close to the correct one)
	int i;
	for (i=0; i<5; i++)
	{
		int idx = (int)(my_random()*pmc_stats.size());
		if (abs(correct_idx-idx)<6)
			continue;

		idxs.push_back(idx);
	}

	sort(idxs.begin(),idxs.end());
	for (i=0; i<idxs.size(); i++)
		if (idxs[i]>=0 && idxs[i]<pmc_stats.size())
			bad_pmc_idxs.push_back(idxs[i]);

}



/*************************************************************************
Tests the performance of precursor mass correction
**************************************************************************/
void PMCSQS_Scorer::test_pmc(Config *config, char *specs_file, int charge, 
							 mass_t min_mass, mass_t max_mass)
{
	BasicSpecReader bsr;
	static QCPeak peaks[5000];

	FileManager fm;
	FileSet fs;
		
	fm.init_from_file(config,specs_file);
	fs.select_files_in_mz_range(fm,min_mass,max_mass,charge);

	const int max_to_read_per_file = 5000;

	const vector<SingleSpectrumFile *>& all_ssf = fs.get_ssf_pointers();
	const int num_samples = (all_ssf.size()<max_to_read_per_file ? all_ssf.size() :
									max_to_read_per_file);
	
	vector<mass_t> org_offsets;
	vector<mass_t> corr_offsets;

	vector<int> ssf_idxs;
	if (num_samples<all_ssf.size())
	{
		choose_k_from_n(num_samples,all_ssf.size(),ssf_idxs);
	}
	else
	{
		int i;
		ssf_idxs.resize(all_ssf.size());
		for (i=0; i<all_ssf.size(); i++)
			ssf_idxs[i]=i;
	}

	vector<SingleSpectrumFile *> ssfs;
	int i;
	for (i=0; i<num_samples; i++)
		ssfs.push_back( all_ssf[ssf_idxs[i]]);
	
	output_pmc_rank_results(fm,charge,ssfs);
}


/***********************************************************************************

Functions for training set.


************************************************************************************/

struct ScanPair {
	ScanPair(int f,int sc, string& se) : file_idx(f), scan(sc), seq(se) {};
	ScanPair(int f,int s) : file_idx(f), scan(s) {};
	ScanPair() : file_idx(-1), scan(-1) {};

	bool operator< (const ScanPair& other) const
	{
		return (file_idx<other.file_idx || 
			    (file_idx == other.file_idx && scan<other.scan));
	}

	bool operator == (const ScanPair& other) const
	{
		return (file_idx == other.file_idx && scan == other.scan);
	}


	int file_idx;
	int scan;
	string seq;
};

void read_idxs_from_file(char *file, vector<ScanPair>& final_pairs, int max_size)
{
	ifstream inp(file,ios::in);
	
	if (! inp.good())
	{	
		cout << "Error opening: " << file << endl;
		exit(1);
	}

	vector<ScanPair> pairs;
	pairs.clear();

	char buff[256];
	while (inp.getline(buff,256))
	{
		istringstream iss(buff);
		int f,s;
		string seq;

		iss >> f >> s >> seq;

		
		if (f>=0 && s>=0)
		{
			if (seq.length()>2)
			{
				pairs.push_back(ScanPair(f,s,seq));
			}
			else
				pairs.push_back(ScanPair(f,s));
		}


	}
	inp.close();

	if (pairs.size() > max_size)
	{
		vector<int> idxs;
		choose_k_from_n(max_size,pairs.size(),idxs);
		final_pairs.resize(max_size);
		int i;
		for (i=0; i<max_size; i++)
			final_pairs[i]=pairs[idxs[i]];
	}
	else
	{
		final_pairs=pairs;
	}

	sort(final_pairs.begin(),final_pairs.end());
}


void create_training_files(Config *config)
{
	char mzxml_list[]={"C:\\Work\\msms5\\PepNovoHQ\\pmcsqs\\HEK293_mzxml_list.txt"};
	char idxs_neg_file[]={"C:\\Work\\msms5\\PepNovoHQ\\pmcsqs\\H40ul_neg_samples.txt"};
//	char idxs1_file[]={"C:\\Work\\msms5\\PepNovoHQ\\pmcsqs\\H40ul_pos_samples.1.txt"};
//	char idxs2_file[]={"C:\\Work\\msms5\\PepNovoHQ\\pmcsqs\\H40ul_pos_samples.2.txt"};
//	char idxs2_file[]={"C:\\Work\\msms5\\PepNovoHQ\\pmcsqs\\Len10_pos_samples.2.txt"};
	char idxs1_file[]={"C:\\Work\\msms5\\PepNovoHQ\\pmcsqs\\sqs_train_pos_samples.1.txt"};
	char idxs2_file[]={"C:\\Work\\msms5\\PepNovoHQ\\pmcsqs\\sqs_train_pos_samples.2.txt"};
	char idxs3_file[]={"C:\\Work\\msms5\\PepNovoHQ\\pmcsqs\\H40ul_pos_samples.3.txt"};

	char out_base[]={"C:\\Work\\msms5\\PepNovoHQ\\pmcsqs\\sqs_train"};
	string out_neg (out_base); 
	string out1=out_neg;
	string out2=out_neg;
	string out3=out_neg;

	out_neg += "_neg.mgf";
	out1 += "_1.mgf";
	out2 += "_2.mgf";
	out3 += "_3.mgf";

	ofstream stream_neg (out_neg.c_str(),ios::out);
	ofstream stream1(out1.c_str(),ios::out);
	ofstream stream2(out2.c_str(),ios::out);
	ofstream stream3(out3.c_str(),ios::out);

	vector<ScanPair> neg_pairs, pairs1,pairs2,pairs3;


	read_idxs_from_file(idxs_neg_file,neg_pairs,12000);
	read_idxs_from_file(idxs1_file,pairs1,12000);
	read_idxs_from_file(idxs2_file,pairs2,12000);
	read_idxs_from_file(idxs3_file,pairs3,8000);

	cout << "Read " << neg_pairs.size() << " neg idxs\n";
	cout << "Read " << pairs1.size() << " pos1 idxs\n";
	cout << "Read " << pairs2.size() << " pos2 idxs\n";
	cout << "Read " << pairs3.size() << " pos3 idxs\n";

	vector<bool> file_inds;
	file_inds.resize(10000,false);
	int i;

	for (i=0; i<neg_pairs.size(); i++)
		file_inds[neg_pairs[i].file_idx]=true;

	for (i=0; i<pairs1.size(); i++)
		file_inds[pairs1[i].file_idx]=true;

	for (i=0; i<pairs2.size(); i++)
		file_inds[pairs2[i].file_idx]=true;

	for (i=0; i<pairs3.size(); i++)
		file_inds[pairs3[i].file_idx]=true;

	
	FileManager fm;
	FileSet fs;

	fm.init_from_list_file(config,mzxml_list,file_inds);
	fs.select_all_files(fm);
	const vector<SingleSpectrumFile *>& all_ssf = fs.get_ssf_pointers();

	

	// read spectra
	BasicSpecReader bsr;
	QCPeak peaks[5000];

	int num_out_neg=0, num_out1=0, num_out2=0, num_out3=0;
	int neg_idx=0,c1_idx=0,c2_idx=0,c3_idx=0;


	for (i=0; i<all_ssf.size(); i++)
	{
		MZXML_single *ssf = (MZXML_single *)all_ssf[i];
		ScanPair ssf_pair(ssf->file_idx,ssf->scan_number);
		string seq="";

		int out_dest=-1;

		while (neg_idx<neg_pairs.size() && neg_pairs[neg_idx]<ssf_pair)
			neg_idx++;

		if (neg_idx<neg_pairs.size() && neg_pairs[neg_idx]==ssf_pair)
			out_dest=0;


		while (c1_idx<pairs1.size() && pairs1[c1_idx]<ssf_pair)
			c1_idx++;
		if (c1_idx<pairs1.size() && pairs1[c1_idx]==ssf_pair)
		{
			seq = pairs1[c1_idx].seq;
			out_dest=1;
		}


		while (c2_idx<pairs2.size() && pairs2[c2_idx]<ssf_pair)
			c2_idx++;
		if (c2_idx<pairs2.size() && pairs2[c2_idx]==ssf_pair)
		{
			seq = pairs2[c2_idx].seq;
			out_dest=2;
		}


		while (c3_idx<pairs3.size() && pairs3[c3_idx]<ssf_pair)
			c3_idx++;
		if (c3_idx<pairs3.size() && pairs3[c3_idx]==ssf_pair)
		{
			seq = pairs3[c3_idx].seq;
			out_dest=3;
		}

		if (out_dest<0)
			continue;

		BasicSpectrum bs;
		bs.num_peaks = bsr.read_basic_spec(config,fm,ssf,peaks);
		bs.peaks = peaks;
		bs.ssf = ssf;

	//	if (out_dest>0)
	//		bs.ssf->peptide.parse_from_string(config,seq);
	
		char name_buff[64];
		if (out_dest==0)
		{
			sprintf(name_buff,"train_neg_%d_%d_%d",num_out_neg,ssf->file_idx,ssf->scan_number);
			bs.ssf->single_name = string(name_buff);
			bs.output_to_mgf(stream_neg,config);
			num_out_neg++;
			continue;
		}

		if (out_dest==1)
		{
			sprintf(name_buff,"train_pos1_%d_%d_%d",num_out1,ssf->file_idx,ssf->scan_number);
			bs.ssf->single_name = string(name_buff);
			bs.output_to_mgf(stream1,config,seq.c_str());
			num_out1++;
			continue;
		}

		if (out_dest==2)
		{
			sprintf(name_buff,"train_pos2_%d_%d_%d",num_out2,ssf->file_idx,ssf->scan_number);
			bs.ssf->single_name = string(name_buff);
			bs.output_to_mgf(stream2,config,seq.c_str());
			num_out2++;
			continue;
		}

		if (out_dest==3)
		{
			sprintf(name_buff,"train_pos3_%d_%d_%d",num_out3,ssf->file_idx,ssf->scan_number);
			bs.ssf->single_name = string(name_buff);
			bs.output_to_mgf(stream3,config,seq.c_str());
			num_out3++;
			continue;
		}
	}

	cout << "Wrote: " << endl;
	cout << "Neg " << num_out_neg << " / " << neg_pairs.size() << endl;
	cout << "Pos1 " << num_out1 << " / " << pairs1.size() << endl;
	cout << "Pos2 " << num_out2 << " / " << pairs2.size() << endl;
	cout << "Pos3 " << num_out3 << " / " << pairs3.size() << endl;

	stream_neg.close();
	stream1.close();
	stream2.close();
	stream3.close();
	
}





void PMCSQS_Scorer::print_spec(const BasicSpectrum& bs) const
{
	cout << bs.ssf->single_name << endl;
	int i;
	for (i=0; i<bs.num_peaks; i++)
	{
		cout << setprecision(2) << fixed << bs.peaks[i].mass << "\t" << bs.peaks[i].intensity << "\t";
		if (curr_spec_iso_levels[i]>0)
			cout << " ISO ";
		if (curr_spec_strong_inds[i])
			cout << " STRONG ";
		cout << endl;
	}
}