quickclusteringspectra.cpp

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

	pkl << charge << endl;


	bool use_scaled=false;
	double scale_factor=1.0;
	double peaks_inten=0;
	if (peaks[0].scaled_intensity>0)
	{
		int i;
		for (i=0; i<peaks.size(); i++)
			peaks_inten += peaks[i].scaled_intensity;

		use_scaled = true;
	}
	else
	{
		int i;
		for (i=0; i<peaks.size(); i++)
			peaks_inten += peaks[i].intensity;
	}

	// avoid very large peak intensity numbers
//	if (total_ion_current > 1000000000)
//		total_ion_current = 1000000000; 

	if (peaks_inten>0 && total_ion_current>0)
	{
		scale_factor = total_ion_current / peaks_inten;
	}

	// filter peaks if there are too many
	if (peaks.size() <= maximum_good_peaks_to_output)
	{	
		int i;
		for (i=0; i<peaks.size(); i++)
		{
			pkl << fixed << setprecision(2) << peaks[i].mass << "\t" << setprecision(1) << 
				scale_factor * (peaks[i].scaled_intensity>0 ? peaks[i].scaled_intensity : peaks[i].intensity);
		//	if (write_peak_count)
		//		pkl << "   " << peaks[i].num_occurences << "/" << peaks[i].max_num_occurences;
			pkl << 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;

			pkl << fixed << setprecision(2) << peaks[i].mass << "\t" << setprecision(1) << 
				scale_factor * (peaks[i].scaled_intensity>0 ? peaks[i].scaled_intensity : peaks[i].intensity);

		//	if (write_peak_count)
		//		pkl << "   " << peaks[i].num_occurences << "/" << peaks[i].max_num_occurences;
			pkl << endl;
		}

	}

	pkl.close();

//	cout << "Wrote: " << file_name << endl;
	
}

void ClusterSpectrum::print_explained_intensity_stats(Peptide& pep) const
{
	int i;
	float total_inten=0;
	vector<mass_t> masses;
	vector<intensity_t> intensities;

	masses.resize(peaks.size());
	intensities.resize(peaks.size());
	for (i=0; i<peaks.size(); i++)
	{
		masses[i]=peaks[i].mass;
		intensities[i]=peaks[i].intensity;
	}

	mass_t pm_with_19 = pep.get_mass()+MASS_OHHH;
	AnnotatedSpectrum as;
	as.read_from_peak_arrays(config,pep,pm_with_19,charge,peaks.size(),
		&masses[0],&intensities[0]);
	as.init_spectrum();

	as.annotate_spectrum(pm_with_19);
	as.print_expected_by();
	
	int len = pep.get_num_aas();
			
	// count b,y
	const int b_frag_idx = config->get_frag_idx_from_label("b");
	const int y_frag_idx = config->get_frag_idx_from_label("y");
	int num_b = as.get_num_observed_frags(b_frag_idx);
	int num_y = as.get_num_observed_frags(y_frag_idx);
	float exp_int = as.get_explianed_intensity();

	cout << "exp: " << exp_int << "  b:" << num_b << "  y:" << num_y <<endl;


}


void BasicSpectrum::print_peaks() const
{
	int i;
	for (i=0; i<this->num_peaks; i++)
		cout << left << setw(5) << i << this->peaks[i].mass << "  " << peaks[i].intensity << endl;
}


// Assumes PTMs are already set
void extractAnnoatedScansFromFiles(Config *config, char *file_list, char *anns, 
								   char *output_file, bool extract_no_process)
{
	vector< vector<int> >    annotation_idxs;
	vector<mzXML_annotation> annotations;
	read_mzXML_annotations(file_list,anns,annotation_idxs,annotations,35000);
	
	cout << "Read annotations: " << annotations.size() << endl;

	FileManager fm;
	fm.init_from_list_file_and_add_annotations(config,file_list,
			annotation_idxs, annotations,true);

	FileSet fs;
	fs.select_all_files(fm,true);
	const vector<SingleSpectrumFile *>& all_ssf = fs.get_ssf_pointers();

	ofstream mgf_stream(output_file,ios::out);
	
	BasicSpecReader bsr;
	QCPeak peaks[5000];
	int i;
	for (i=0; i<all_ssf.size(); i++)
	{
		BasicSpectrum bs;
		MZXML_single *ssf = (MZXML_single *)all_ssf[i];
	
		bs.peaks = peaks;
		bs.ssf = ssf;
		bs.ssf->charge   = ssf->charge;

		bs.num_peaks = bsr.read_basic_spec(config,fm,ssf,peaks,extract_no_process);
		
		if (ssf->scan_number<0)
		{
			cout << "Error: no scan number read from MGF!!!" << endl;
			exit(1);
		}

		bs.output_to_mgf(mgf_stream,config);
	}
	mgf_stream.close();
}





void read_mzXML_annotations_to_map(char *ann_file, 
								   map<mzXML_annotation,int>& ann_map)
{

	int i;
	char buff[256];


	FILE *ann_stream = fopen(ann_file,"r");
	if (! ann_stream)
	{
		cout << "Error: couldn't open annotation files for run!: " << ann_file << endl;
		exit(1);
	}



	int anns=0;
	i=0;
	while (fgets(buff,256,ann_stream))
	{
		int file_idx=-1, mzXML_idx=-1; 
		char only_peptide[128];
		int scan=-1,charge=0;

		if (sscanf(buff,"%d %d %d %d %s",&file_idx,&mzXML_idx,&scan,&charge,only_peptide)<5)
			continue;

		mzXML_annotation ann;
		ann.charge =charge;
		ann.scan = scan;
		ann.mzXML_file_idx = file_idx;
		ann.pep = only_peptide;


		ann_map.insert(make_pair(ann,1));
		anns++;
	}
	cout << "Read " << anns << " annotations..." << endl;
}


// reads an annotation file. mzXML
// saves the annotation for each file number and scan number
// assume max scan number is 30000
// file_idx  mzXML_file_idx  scan peptide
// 169 01 8854 2 VAQGVSGAVQDK
void read_mzXML_annotations(char *mzXML_list, 
							char *ann_file, 
							vector< vector<int> >& annotation_idxs, 
							vector<mzXML_annotation>& annotations,
							int max_ann_size) 
{
	int i;
	char buff[256];
	FILE *mzxml_stream = fopen(mzXML_list,"r");
	if (! mzxml_stream)
	{
		cout << "Error: couldn't open annotation file for mzXML run!: " << mzXML_list << endl;
		exit(1);
	}

	int count=0;
	while (fgets(buff,256,mzxml_stream))
	{
		count++;
	}
	fclose(mzxml_stream);

	annotation_idxs.resize(count+1);
	for (i=0; i<=count; i++)
		annotation_idxs[i].resize(max_ann_size,-1);

	FILE *ann_stream = fopen(ann_file,"r");
	if (! ann_stream)
	{
		cout << "Error: couldn't open annotation files for run!: " << ann_file << endl;
		exit(1);
	}

/*	b-total-try-2nd-digest-b-400ug-2D34-121505-LTQ2-19.mzXML'...
20989 spectra...
255 Parse spectra from 'C:/Work/Data/Briggs\H293b-total-try-2nd-digest-b-400ug-2D34-121505-LTQ2\H293
b-total-try-2nd-digest-b-400ug-2D34-121505-LTQ2-20.mzXML'...
20712 spectra...
256 Parse spectra from 'C:/Work/Data/Briggs\H293b-total-try-2nd-digest-b-400ug-2D34-121505-LTQ2\H293
b-total-try-2nd-digest-b-400ug-2D34-121505-LTQ2-21.mzXML'...
20603 spectra...*/

	

	i=0;
	while (fgets(buff,256,ann_stream))
	{
		int file_idx=-1, mzXML_idx=-1; 
		char only_peptide[128];
		int scan=-1,charge=0;

		if (sscanf(buff,"%d %d %d %d %s",&file_idx,&mzXML_idx,&scan,&charge,only_peptide)<5)
			continue;

		mzXML_annotation ann;
		ann.charge =charge;
		ann.mzXML_file_idx = file_idx;
		ann.pep = only_peptide;


		annotation_idxs[file_idx][scan]=annotations.size();
		annotations.push_back(ann);
		
	//	cout << file_idx << " : >> " << scan << "   " << only_peptide << " " << charge << endl;

	//	if (++i>=200)
	//		break;
	
	}
}






void FileManager::init_from_dat_list_extract_only_annotated(Config *config, 
						char* dat_list_file, char *ann_file)
{

	int i;

	// read annotations
	map<mzXML_annotation,int> ann_map;
	read_mzXML_annotations_to_map(ann_file,ann_map);

	vector<string> list;
	read_paths_into_list(dat_list_file,list);
	dat_files.clear();
	for (i=0; i<list.size(); i++)
    {
		if (list[i][0] == '#')
			continue;

		DAT_file dat;
		dat.dat_name =list[i];

		dat.initial_read(config,dat_files.size());

		cout << dat.dat_name << " .. ";

		// change the single spectrum pointers in the mgf file record
		// to include only those that have a mass that is in the permitted range

		vector<DAT_single> good_singles;
		int j;
		for (j=0; j<dat.single_spectra.size(); j++)
		{
			int mzxml_file_idx = dat.single_spectra[j].mzxml_file_idx;
			int scan_number = dat.single_spectra[j].scan_number;

			map<mzXML_annotation,int>::const_iterator it;
			mzXML_annotation ann_pos;

			ann_pos.mzXML_file_idx = mzxml_file_idx;
			ann_pos.scan = scan_number;

			it = ann_map.find(ann_pos);

			if (it == ann_map.end())
				continue;

			Peptide pep;
			int charge = it->first.charge;
			const string& pep_str = it->first.pep;

			pep.parse_from_string(config,pep_str);
			pep.calc_mass(config);
			mass_t m_over_z = (pep.get_mass()+19.0183 + (mass_t)charge)/charge;

			if (fabs(m_over_z - dat.single_spectra[j].m_over_z)>7.0)
			{
				cout << "Error: mismatch between ann " << pep_str << " " << m_over_z << endl;
				cout << "       and dat  " << mzxml_file_idx  << ", " << scan_number << "  " << dat.single_spectra[j].m_over_z << endl;

				cout << "Mass Cys: " << config->get_aa2mass()[Cys] << endl;
				continue;
			}
					
			dat.single_spectra[j].peptide.parse_from_string(config,
					it->first.pep);
		
			good_singles.push_back(dat.single_spectra[j]);
		}

		dat.single_spectra = good_singles;
		dat_files.push_back(dat);

		cout << good_singles.size() << " ..." << endl;

	}

}



void extract_annotate_scans_from_dat(Config *config, char *dat_list, char *anns_file, 
									 char *out_name)
{
	FileManager fm;
	fm.init_from_dat_list_extract_only_annotated(config,dat_list,anns_file);

	FileSet fs;
	fs.select_all_files(fm,true);
	const vector<SingleSpectrumFile *>& all_ssf = fs.get_ssf_pointers();

	cout << "Processing " << all_ssf.size() << " spectra..." << endl;


	vector<FILE *> len_streams,size_streams;
	vector<int> len_counts, size_counts;

	const int max_len = 60;
	const int max_size_idx = 65;

	len_streams.resize(max_len+1,NULL);
	size_streams.resize(max_size_idx+1,NULL);

	len_counts.resize(max_len+1,0);
	size_counts.resize(max_size_idx+1,0);

	
	BasicSpecReader bsr;
	QCPeak peaks[5000];
	int i;
	for (i=0; i<all_ssf.size(); i++)
	{
		BasicSpectrum bs;
		DAT_single *ssf = (DAT_single *)all_ssf[i];
	
		bs.peaks = peaks;
		bs.ssf = ssf;
		bs.ssf->charge   = ssf->charge;

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

		const int len_idx = bs.ssf->peptide.get_num_aas();

		bs.ssf->peptide.calc_mass(config);
		const int size_idx = (int)((bs.ssf->peptide.get_mass()+MASS_OHHH)/100);

		if (len_idx>max_len || size_idx>max_size_idx)
			continue;
		
		if (ssf->scan_number<0)
		{
			cout << "Error: no scan number read from MGF!!!" << endl;
			exit(1);
		}

		char title[32];

		sprintf(title,"%d_%d",ssf->mzxml_file_idx,ssf->scan_number);
		ssf->single_name = string(title);



		ofstream len_stream,size_stream;

		if (! len_streams[len_idx])
		{
			char len_name[256];
			sprintf(len_name,"%s_%d.mgf",out_name,len_idx);
			len_streams[len_idx] = fopen(len_name,"w");
			cout << "open: " << len_name << endl;
		}
	

		if (! size_streams[size_idx])
		{
			char size_name[256];
			sprintf(size_name,"%s_%d.mgf",out_name,size_idx*100);
		
			size_streams[size_idx] = fopen(size_name,"w");

			cout << "open: " << size_name << endl;
		}
	
		bs.output_to_mgf(len_streams[len_idx],config);
		bs.output_to_mgf(size_streams[size_idx],config);

		len_counts[len_idx]++;
		size_counts[size_idx]++;
	}


	for (i=0; i<len_streams.size(); i++)
		if (len_streams[i])
			fclose(len_streams[i]);

	for (i=0; i<size_streams.size(); i++)
		if (size_streams[i])
			fclose(size_streams[i]);

	char len_sum_name[256],size_sum_name[256];
	sprintf(len_sum_name,"%s_len_sum.txt",out_name);
	sprintf(size_sum_name,"%s_size_sum.txt",out_name);

	ofstream len_sum(len_sum_name,ios::out);
	for (i=0; i<=max_len; i++)
		if (len_counts[i]>0)
			len_sum << i << "\t" << len_counts[i] << endl;
	len_sum.close();

	ofstream size_sum(size_sum_name,ios::out);
	for (i=0; i<=max_size_idx; i++)
		if (size_counts[i]>0)
			size_sum << i*100 << "\t" << size_counts[i] << endl;
	size_sum.close();
}



void convert_dat_to_mgf(Config *config, 
						char *dat_list, 
						char *out_name, 
						char *out_dir,
						char *anns_file)
{
	QCOutputter qco_all,qco_labeled;

	qco_all.init(string(out_name) + "_all",out_dir);
	qco_labeled.init(string(out_name) + "_labeled",out_dir);

	map<mzXML_annotation,int> ann_map;
	bool use_map = false;
	if (anns_file)
	{
		read_mzXML_annotations_to_map(anns_file,ann_map);
		if (ann_map.size()>1)
			use_map=true;
	}

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

	BasicSpecReader bsr;
	QCPeak peaks[5000];

	cout << "Processing " << all_ssf.size() << " spectra..." << endl;
	int num_wrote = 0;
	int num_annotated = 0;
	int bad_anns= 0;
	int i;
	for (i=0; i<all_ssf.size(); i++)
    {
		DAT_single *ssf = (DAT_single *)all_ssf[i];
		const int mzxml_file_idx = ssf->mzxml_file_idx;
		const int scan_number = ssf->scan_number;

		map<mzXML_annotation,int>::const_iterator it;
		mzXML_annotation ann_pos;

		ann_pos.mzXML_file_idx = mzxml_file_idx;
		ann_pos.scan = scan_number;

		it = ann_map.find(ann_pos);

		bool has_ann = false;
		if (it != ann_map.end())
		{
			ssf->peptide.parse_from_string(config,it->first.pep);
			ssf->peptide.calc_mass(config);
			const int charge = it->first.charge;
			mass_t exp_m_over_z = (ssf->peptide.get_mass() + 18.0 + charge)/charge;
			if (fabs(ssf->m_over_z-exp_m_over_z)>6.0)
			{
				bad_anns++;
				has_ann=false;
				ssf->peptide.clear();
				continue;
			}
			{
				num_annotated++;
				has_ann=true;
			}
		}
	
		BasicSpectrum bs;
		
		bs.peaks = peaks;
		bs.ssf = ssf;
		bs.ssf->charge   = ssf->charge;

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

		if (bs.num_peaks<10)
			continue;

		char title[32];

		sprintf(title,"%d_%d",ssf->mzxml_file_idx,ssf->scan_number);
		ssf->single_name = string(title);

		qco_all.output_basic_spectrum_to_mgf(bs,config);
		if (it != ann_map.end())
			qco_labeled.output_basic_spectrum_to_mgf(bs,config);
		num_wrote++;
	}

	cout << "Wrote " << num_wrote << " spectra to MGF" << endl;
	cout << num_annotated << " were annotated. " << endl;
	cout << "Found " << bad_anns << " bad anns..." << endl;
}