pmcsqs.cpp

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

		}

		// find indicators for min tolerance and max pairs
		float tol_pairs =           POS_INF;
		float tol_strong_pairs =    POS_INF;
		float tol_c2_pairs =        POS_INF;
		float tol_c2_strong_pairs = POS_INF;

		int   idx_pairs=0;
		int	  idx_strong_pairs=0;
		int   idx_c2_pairs=0;
		int   idx_c2_strong_pairs=0;

		for (i=0; i<curr_spec_rank_pmc_tables[charge].size(); i++)
		{
			const PMCRankStats& curr_stats = curr_spec_rank_pmc_tables[charge][i];

			if (curr_stats.num_frag_pairs == maximal.num_frag_pairs &&
				curr_stats.mean_offset_pairs < tol_pairs)
			{
				idx_pairs=i;
				tol_pairs=curr_stats.mean_offset_pairs;
			}

			if (curr_stats.num_strong_frag_pairs == maximal.num_strong_frag_pairs &&
				curr_stats.mean_offset_c2_strong_pairs < tol_strong_pairs)
			{
				idx_strong_pairs=i;
				tol_strong_pairs=curr_stats.mean_offset_c2_strong_pairs;
			}

			if (curr_stats.num_c2_frag_pairs == maximal.num_c2_frag_pairs &&
				curr_stats.mean_offset_c2_pairs < tol_c2_pairs)
			{
				idx_c2_pairs=i;
				tol_c2_pairs=curr_stats.mean_offset_c2_pairs;
			}

			if (curr_stats.num_strong_c2_frag_pairs == maximal.num_strong_c2_frag_pairs &&
				curr_stats.mean_offset_c2_strong_pairs < tol_c2_strong_pairs)
			{
				idx_c2_strong_pairs=i;
				tol_c2_strong_pairs=curr_stats.mean_offset_c2_strong_pairs;
			}

		}

		curr_spec_rank_pmc_tables[charge][idx_pairs].ind_pairs_with_min_tol=true;
		curr_spec_rank_pmc_tables[charge][idx_strong_pairs].ind_strong_pairs_with_min_tol=true;
		curr_spec_rank_pmc_tables[charge][idx_c2_pairs].ind_c2_pairs_with_min_tol=true;
		curr_spec_rank_pmc_tables[charge][idx_c2_strong_pairs].ind_c2_strong_pairs_with_min_tol=true;


		static vector<float> log_distances;
		if (log_distances.size()<curr_spec_rank_pmc_tables[charge].size())
		{
			log_distances.resize(curr_spec_rank_pmc_tables[charge].size(),0);
			int i;
			for (i=1; i<log_distances.size(); i++)
				log_distances[i]=log(1.0+(float)i);
		}

		for (i=0; i<curr_spec_rank_pmc_tables[charge].size(); i++)
		{
			PMCRankStats& curr_stats = curr_spec_rank_pmc_tables[charge][i];
			curr_stats.log_dis_from_pairs_min_tol = log_distances[abs(i-idx_pairs)];
			curr_stats.log_dis_from_strong_pairs_min_tol = log_distances[abs(i-idx_strong_pairs)];
			curr_stats.log_dis_from_c2_pairs_min_tol = log_distances[abs(i-idx_c2_pairs)];
			curr_stats.log_dis_from_c2_strong_pairs_min_tol = log_distances[abs(i-idx_c2_strong_pairs)];
		}

	}
}
	


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

******************************************************************************************/
void PMCSQS_Scorer::get_sqs_features_from_pmc_tables(const BasicSpectrum& bs,
						vector< vector<float> >& sqs_features) const
{
	const mass_t org_m_over_z = bs.ssf->m_over_z;
	float max_num_strong_pairs=0;
	int   best_table_idx=0;
	mass_t mz_diff = 99999;

	if (sqs_features.size() != max_model_charge +1)
		sqs_features.resize(max_model_charge+1);

	int charge;
	for (charge=1; charge<=max_model_charge; charge++)
	{
		// find entry which has the maximal number of strong pairs, while maining the minimial
		// m/z distance shift from the original
		int i;
		for (i=1; i<this->curr_spec_rank_pmc_tables[charge].size(); i++)
		{
			const float curr_num_strong = curr_spec_rank_pmc_tables[charge][i].num_strong_frag_pairs;
			
			if (curr_num_strong>=max_num_strong_pairs)
			{
				const mass_t distance = fabs(curr_spec_rank_pmc_tables[charge][i].m_over_z - org_m_over_z);
				if (curr_num_strong == max_num_strong_pairs && distance>= mz_diff)
					continue;
				
				max_num_strong_pairs = curr_num_strong;
				mz_diff = distance;
				best_table_idx = i;
			}
		}

		const PMCRankStats& best_stats = curr_spec_rank_pmc_tables[charge][best_table_idx];

		sqs_features[charge].clear();

		sqs_features[charge].push_back(best_stats.num_frag_pairs);

		sqs_features[charge].push_back(best_stats.num_strong_frag_pairs);

		sqs_features[charge].push_back(best_stats.num_c2_frag_pairs);

		sqs_features[charge].push_back(best_stats.num_strong_c2_frag_pairs);
	}
}



void PMCSQS_Scorer::compute_sqs_cum_stats_for_ided(Config *config, char *list)
{
	const float max_prob=0.5;
	const int max_bin_idx = int(max_prob*100)+1;
	vector< vector< vector<int> > > bin_counts; //charge / size_idx / bin
	vector< vector< int > > totals;

	bin_counts.resize(4);
	totals.resize(4);
	int i;
	for (i=0; i<4; i++)
	{
		int n = this->sqs_mass_thresholds.size() + 1;
		totals[i].resize(n,0);
		bin_counts[i].resize(n);
		int j;
		for (j=0; j<n; j++)
			bin_counts[i][j].resize(max_bin_idx+1,0);
	}

	// prefrom sqs on all

		

	FileManager fm;
	FileSet fs;

	fm.init_from_list_file(config,list);
	fs.select_all_files(fm);

	const int max_to_read_per_file = 100000;

	const vector<SingleSpectrumFile *>& all_ssf = fs.get_ssf_pointers();

	BasicSpecReader bsr;
	static QCPeak peaks[5000];

	int low_count = 0;
	for (i=0; i<all_ssf.size(); i++)
	{
		SingleSpectrumFile* ssf = all_ssf[i];
		BasicSpectrum bs;
	
		bs.num_peaks = bsr.read_basic_spec(config,fm,ssf,peaks);
		bs.peaks = peaks;
		bs.ssf = ssf;

		const int org_charge=ssf->peptide.calc_charge(ssf->m_over_z);
		const int size_idx = get_sqs_size_idx(ssf->m_over_z);
		int charge=0;
		float prob=get_sqs_for_spectrum(config,bs,&charge);

		int bin_idx=int(100*prob);
		if (bin_idx>max_bin_idx)
			bin_idx = max_bin_idx;

		totals[org_charge][size_idx]++;
		int j;
		for (j=0; j<=bin_idx; j++)
			bin_counts[org_charge][size_idx][j]++;

		if (prob<0.1)
		{
			cout << ++low_count << "\t" << ssf->single_name << "\t" << org_charge << " : " << 
				charge << " --> " << prob << endl;
		}
	}

	int c;
	for (c=0; c<4; c++)
	{
		int specs=0;
		int i;
		for (i=0; i<totals[c].size(); i++)
			specs+=totals[c][i];
		if (specs<100)
			continue;

		cout << endl << "CHARGE " << c << endl;
		for (i=0; i<totals[c].size(); i++)
		{
			if (totals[c][i]<50)
				continue;
			cout << "SIZE " << i << " ( charge " << c << " )" << endl;
			int j;
			cout << setprecision(4);
			for (j=0; j<=max_bin_idx; j++)
			{
				cout << j*0.01 << "\t" << bin_counts[c][i][j] << "\t" << 
					bin_counts[c][i][j]/(float)totals[c][i] << endl;
			}
			cout << endl;
		}
	}
}






void PMCSQS_Scorer::create_filtered_peak_list_for_sqs(
									  QCPeak *org_peaks, int num_org_peaks,
									  QCPeak *new_peaks, int& num_new_peaks) const
{
	static vector<MassInten> peak_list;
	static int peak_list_size =0;
	int i;

	if (num_org_peaks>peak_list_size)
	{
		peak_list_size = (int)(num_org_peaks * 1.5);
		if (peak_list_size<2000)
			peak_list_size = 2000;

		peak_list.resize(peak_list_size);
	}

	// copy org_peaks to the temporary peak_list
	int f_idx=0;
	for (i=0; i<num_org_peaks; i++)
	{
		peak_list[i].mass=org_peaks[i].mass;
		peak_list[i].intensity=org_peaks[i].intensity;
	}

	const mass_t tolerance = config->get_tolerance();
	const mass_t join_tolerance = (tolerance < 0.05 ? tolerance : 0.5 * tolerance);
	int p_idx=0;
	i=1;
	while (i<num_org_peaks)
	{
		if (peak_list[i].mass - peak_list[p_idx].mass<=join_tolerance)
		{
			intensity_t inten_sum = peak_list[i].intensity + peak_list[p_idx].intensity;
			mass_t new_mass = (peak_list[i].intensity * peak_list[i].mass + 
							   peak_list[p_idx].intensity * peak_list[p_idx].mass ) / inten_sum;

			peak_list[p_idx].mass = new_mass;
			peak_list[p_idx].intensity = inten_sum;	
		}
		else
		{
			peak_list[++p_idx]=peak_list[i];
		}
		i++;
	}
	int num_peaks = p_idx+1;


	// filter low intensity noise
	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();
	const int max_peak_idx = num_peaks -1;
	int min_idx=1;
	int max_idx=1;
	p_idx =1;

	
	new_peaks[0].mass      = peak_list[0].mass;
	new_peaks[0].intensity = peak_list[0].intensity;
	f_idx=1;

	// 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)
		{
			new_peaks[f_idx].mass = peak_list[i].mass;
			new_peaks[f_idx].intensity = peak_list[i].intensity;
			f_idx++;
			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)
		{
			new_peaks[f_idx].mass = peak_list[i].mass;
			new_peaks[f_idx].intensity = peak_list[i].intensity;
			f_idx++;
		}
	}
	new_peaks[f_idx].mass = peak_list[i].mass;
	new_peaks[f_idx].intensity = peak_list[i].intensity;
	f_idx++;


	num_new_peaks = f_idx;

	// normalize intensities

	if (1)
	{
		intensity_t total_inten=0;

		for (i=1; i<num_new_peaks; i++)
			total_inten+=new_peaks[i].intensity;

		const mass_t one_over_total_inten = (1000.0 / total_inten);

		for (i=1; i<num_new_peaks; i++)
			new_peaks[i].intensity *= one_over_total_inten; 
	}
}



/*****************************************************************
Creates for each input file an mgf file that holds the spectra
that passed quality filtering does not correct PM and charge in the 
mgf files.
******************************************************************/
void PMCSQS_Scorer::output_filtered_spectra_to_mgfs(
									 Config *config,
									 const vector<string>& files,
									 char *out_dir,
									 float filter_prob, 
									 int& total_num_written, 
									 int& total_num_read)
{
	total_num_read = 0;
	total_num_read = 0;
	int f;
	for (f=0; f<files.size(); f++)
	{
		const char *spectra_file = files[f].c_str();
		FileManager fm;
		FileSet fs;
		BasicSpecReader bsr;

		string fname, mgf_name, map_name;
		get_file_name_without_extension(files[f],fname);

		mgf_name = string(out_dir) + "/" + fname + "_fil.mgf";
		map_name = string(out_dir) + "/" + fname + "_map.txt";

		///////////////////////////////////////////////
		// Quick read, get all pointers to begining of spectra
		if (get_file_extension_type(files[f]) != MZXML)
		{
			fm.init_from_file(config,spectra_file);
		}
		else
			fm.init_and_read_single_mzXML(config,spectra_file,f);

		fs.select_all_files(fm);

		const vector<SingleSpectrumFile *>& all_ssf = fs.get_ssf_pointers();
		int sc;
		int  num_spec_written=0;
		bool first=true;
		ofstream out_stream, map_stream;
		for (sc=0; sc<all_ssf.size(); sc++)
		{
			static vector<QCPeak> peaks;
			SingleSpectrumFile *ssf = all_ssf[sc];
			
			if (peaks.size()<ssf->num_peaks)
			{
				int new_size = ssf->num_peaks*2;
				if (new_size<2500)
					new_size=2500;
				peaks.resize(new_size);
			}

			// read without processing peaks
			const int num_peaks = bsr.read_basic_spec(config,fm,ssf,&peaks[0],false,true);
			ssf->file_idx = f;

			BasicSpectrum bs;
			bs.peaks = &peaks[0];
			bs.num_peaks = num_peaks;
			bs.ssf = all_ssf[sc];

			int max_charge;
			float prob = get_sqs_for_spectrum(config,bs,&max_charge);
			if (prob<filter_prob)
			{
				continue;
			}

			if (first)
			{
				out_stream.open(mgf_name.c_str(),ios::out);
				map_stream.open(map_name.c_str(),ios::out);

				cout << "Filtering spectra to minumum quality score: " << filter_prob << endl;
				cout << "Writing spectra info to:" << endl;
				cout << mgf_name << endl << map_name << endl;

				if (! out_stream.is_open() || ! out_stream.good())
				{
					cout << "Error: couldn\'t open for out mgf stream for writing: " <<
						endl << mgf_name << endl;
					exit(1);
				}
				first = false;
			}

			char single_name[64];
			sprintf(single_name,"%d:%d",f,bs.ssf->get_scan());
			bs.ssf->single_name = single_name;
			bs.output_to_mgf(out_stream,config);
			if (prob>1.0)
				prob=1.0;
			map_stream << num_spec_written++ << "\t" << all_ssf[sc]->get_scan() << "\t" << fixed << prob << endl;

			
		}
		out_stream.close();
		map_stream.close();

		total_num_read+= all_ssf.size();
		total_num_written += num_spec_written;
	}
}