ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/gclib/tophat_cpp/juncs_db.cpp

#	Line 26 | Line 26
26		#include <seqan/sequence.h>
27		#include <seqan/find.h>
28		#include <seqan/file.h>
29	+	#include <seqan/modifier.h>
30		#include <getopt.h>
31
32		#include "common.h"
#	Line 34 | Line 35
35		#include "junctions.h"
36		#include "insertions.h"
37		#include "deletions.h"
38	+	#include "fusions.h"
39	+
40
41		using namespace std;
42		using namespace seqan;
#	Line 98 | Line 101
101
102		template<typename TStr>
103		void print_splice(const Junction& junction,
104	<	int read_len,
105	<	const string& tag,
106	<	TStr& ref_str,
107	<	const string& ref_name,
108	<	ostream& splice_db)
104	>	int read_len,
105	>	const string& tag,
106	>	TStr& ref_str,
107	>	const string& ref_name,
108	>	ostream& splice_db)
109		{
110		// daehwan - this is tentative, let's think about this more :)
111		// int half_splice_len = read_len - min_anchor_len;
112		int half_splice_len = read_len;
110	–
113		size_t left_start, right_start;
114		size_t left_end, right_end;
115
#	Line 134 | Line 136
136		}
137		}
138
139	+	template<typename TStr>
140	+	void print_fusion(const Fusion& fusion,
141	+	int read_len,
142	+	TStr& left_ref_str,
143	+	TStr& right_ref_str,
144	+	const char* left_ref_name,
145	+	const char* right_ref_name,
146	+	ostream& fusion_db)
147	+	{
148	+	int half_splice_len = read_len - min_anchor_len;
149	+
150	+	size_t left_start, right_start;
151	+	size_t left_end, right_end;
152	+
153	+	if (fusion.left >= 0 && fusion.left < length(left_ref_str) &&
154	+	fusion.right >= 0 && fusion.right < length(right_ref_str))
155	+	{
156	+	if (fusion.dir == FUSION_FF || fusion.dir == FUSION_FR)
157	+	{
158	+	left_start = fusion.left + 1 >= half_splice_len ? fusion.left - half_splice_len + 1 : 0;
159	+	left_end = left_start + half_splice_len;
160	+	}
161	+	else
162	+	{
163	+	left_start = fusion.left;
164	+	left_end = left_start + half_splice_len < length(left_ref_str) ? left_start + half_splice_len : length(left_ref_str);
165	+	}
166	+
167	+	if (fusion.dir == FUSION_FF || fusion.dir == FUSION_RF)
168	+	{
169	+	right_start = fusion.right;
170	+	right_end = right_start + half_splice_len < length(right_ref_str) ? right_start + half_splice_len : length(right_ref_str);
171	+	}
172	+	else
173	+	{
174	+	right_end = fusion.right + 1;
175	+	right_start = right_end >= half_splice_len ? right_end - half_splice_len : 0;
176	+	}
177	+
178	+	seqan::Dna5String left_splice = infix(left_ref_str, left_start, left_end);
179	+	seqan::Dna5String right_splice = infix(right_ref_str, right_start, right_end);
180	+
181	+	if (fusion.dir == FUSION_RF || fusion.dir == FUSION_RR)
182	+	{
183	+	seqan::convertInPlace(left_splice, seqan::FunctorComplement<Dna>());
184	+	seqan::reverseInPlace(left_splice);
185	+
186	+	left_start = left_end - 1;
187	+	}
188	+
189	+	if (fusion.dir == FUSION_FR || fusion.dir == FUSION_RR)
190	+	{
191	+	seqan::convertInPlace(right_splice, seqan::FunctorComplement<Dna>());
192	+	seqan::reverseInPlace(right_splice);
193	+
194	+	right_end = right_start - 1;
195	+	}
196	+
197	+	const char* dir = "ff";
198	+	if (fusion.dir == FUSION_FR)
199	+	dir = "fr";
200	+	else if (fusion.dir == FUSION_RF)
201	+	dir = "rf";
202	+	else if (fusion.dir == FUSION_RR)
203	+	dir = "rr";
204	+
205	+	fusion_db << ">" << left_ref_name << "-" << right_ref_name << "|"
206	+	<< left_start << "|"
207	+	<< fusion.left << "-" << fusion.right << "|"
208	+	<< right_end << "|fus|" << dir <<  endl;
209	+
210	+	fusion_db << left_splice << right_splice << endl;
211	+	}
212	+	}
213	+
214
215		/**
216		* Parse an int out of optarg and enforce that it be at least 'lower';
#	Line 179 | Line 256
256		//      return 0;
257		//}
258
259	+	void get_seqs(istream& ref_stream,
260	+	RefSequenceTable& rt,
261	+	bool keep_seqs = true)
262	+	{
263	+	while(ref_stream.good() &&
264	+	!ref_stream.eof())
265	+	{
266	+	RefSequenceTable::Sequence* ref_str = new RefSequenceTable::Sequence();
267	+	string name;
268	+	readMeta(ref_stream, name, Fasta());
269	+	string::size_type space_pos = name.find_first_of(" \t\r");
270	+	if (space_pos != string::npos)
271	+	{
272	+	name.resize(space_pos);
273	+	}
274	+	fprintf(stderr, "\tLoading %s...", name.c_str());
275	+	seqan::read(ref_stream, *ref_str, Fasta());
276	+	fprintf(stderr, "done\n");
277	+	rt.get_id(name, keep_seqs ? ref_str : NULL, 0);
278	+	if (!keep_seqs)
279	+	delete ref_str;
280	+	}
281	+	}
282	+
283		void driver(const vector<FILE*>& splice_coords_files,
284	<	const vector<FILE*>& insertion_coords_files,
285	<	const vector<FILE*>& deletion_coords_files,
286	<	ifstream& ref_stream)
284	>	const vector<FILE*>& insertion_coords_files,
285	>	const vector<FILE*>& deletion_coords_files,
286	>	const vector<FILE*>& fusion_coords_files,
287	>	ifstream& ref_stream)
288		{
289		char splice_buf[2048];
290	<	RefSequenceTable rt(true);
290	>	RefSequenceTable rt(sam_header, true);
291	>	get_seqs(ref_stream, rt, true);
292	>
293		JunctionSet junctions;
294		for (size_t i = 0; i < splice_coords_files.size(); ++i)
295		{
#	Line 316 | Line 420
420		}
421		}
422
423	+	std::set<Fusion> fusions;
424	+	for(size_t i=0; i < fusion_coords_files.size(); ++i){
425	+	FILE* fusion_coords = fusion_coords_files[i];
426	+	if(!fusion_coords){
427	+	continue;
428	+	}
429	+	while(fgets(splice_buf, 2048, fusion_coords)){
430	+	char* nl = strrchr(splice_buf, '\n');
431	+	char* buf = splice_buf;
432	+	if (nl) *nl = 0;
433	+
434	+	char* ref_name1 = strsep((char**)&buf, "\t");
435	+	char* scan_left_coord = strsep((char**)&buf, "\t");
436	+	char* ref_name2 = strsep((char**)&buf, "\t");
437	+	char* scan_right_coord = strsep((char**)&buf, "\t");
438	+	char* scan_dir = strsep((char**)&buf, "\t");
439
440	<	typedef RefSequenceTable::Sequence Reference;
441	<
442	<	while(ref_stream.good() &&
443	<	!ref_stream.eof())
444	<	{
445	<	Reference ref_str;
446	<	string name;
447	<
448	<	readMeta(ref_stream, name, Fasta());
449	<	string::size_type space_pos = name.find_first_of(" \t\r");
450	<	if (space_pos != string::npos)
451	<	{
452	<	name.resize(space_pos);
453	<	}
454	<
455	<	read(ref_stream, ref_str, Fasta());
456	<
457	<	uint32_t refid = rt.get_id(name, NULL, 0);
458	<	Junction dummy_left(refid, 0, 0, true);
339	<	Junction dummy_right(refid, VMAXINT32, VMAXINT32, true);
340	<
341	<	pair<JunctionSet::iterator, JunctionSet::iterator> r;
342	<	r.first = junctions.lower_bound(dummy_left);
343	<	r.second = junctions.upper_bound(dummy_right);
344	<
345	<	JunctionSet::iterator itr = r.first;
346	<
347	<	while(itr != r.second && itr != junctions.end())
348	<	{
349	<	print_splice(itr->first, read_length, itr->first.antisense ? "GTAG|rev" : "GTAG|fwd", ref_str, name, cout);
350	<	++itr;
440	>	if (!ref_name1 || !scan_left_coord || !ref_name2 || !scan_right_coord || !scan_dir)
441	>	{
442	>	fprintf(stderr,"Error: malformed insertion coordinate record\n");
443	>	exit(1);
444	>	}
445	>
446	>	uint32_t ref_id1 = rt.get_id(ref_name1, NULL, 0);
447	>	uint32_t ref_id2 = rt.get_id(ref_name2, NULL, 0);
448	>	uint32_t left_coord = atoi(scan_left_coord);
449	>	uint32_t right_coord = atoi(scan_right_coord);
450	>	uint32_t dir = FUSION_FF;
451	>	if (strcmp(scan_dir, "fr") == 0)
452	>	dir = FUSION_FR;
453	>	else if(strcmp(scan_dir, "rf") == 0)
454	>	dir = FUSION_RF;
455	>	else if(strcmp(scan_dir, "rr") == 0)
456	>	dir = FUSION_RR;
457	>
458	>	fusions.insert(Fusion(ref_id1, ref_id2, left_coord, right_coord, dir));
459		}
460		}
461
354	–
355	–	ref_stream.clear();
356	–	ref_stream.seekg(0, ios::beg);
357	–
358	–
359	–	while(ref_stream.good() &&
360	–	!ref_stream.eof())
462		{
463	<	Reference ref_str;
464	<	string name;
465	<
466	<	readMeta(ref_stream, name, Fasta());
467	<	string::size_type space_pos = name.find_first_of(" \t\r");
468	<	if (space_pos != string::npos)
469	<	{
470	<	name.resize(space_pos);
370	<	}
371	<
372	<	read(ref_stream, ref_str, Fasta());
373	<
374	<	uint32_t refid = rt.get_id(name, NULL,0);
375	<	Deletion dummy_left(refid, 0, 0, true);
376	<	Deletion dummy_right(refid, VMAXINT32, VMAXINT32, true);
377	<
378	<	pair<std::set<Deletion>::iterator, std::set<Deletion>::iterator> r;
379	<	r.first = deletions.lower_bound(dummy_left);
380	<	r.second = deletions.upper_bound(dummy_right);
381	<
382	<	std::set<Deletion>::iterator itr = r.first;
383	<
384	<	while(itr != r.second && itr != deletions.end())
385	<	{
386	<	print_splice((Junction)*itr, read_length, itr->antisense ? "del|rev" : "del|fwd", ref_str, name, cout);
387	<	++itr;
388	<	}
463	>	JunctionSet::iterator itr = junctions.begin();
464	>	while(itr != junctions.end())
465	>	{
466	>	RefSequenceTable::Sequence* ref_str = rt.get_seq(itr->first.refid);
467	>	const char* name = rt.get_name(itr->first.refid);
468	>	print_splice(itr->first, read_length, itr->first.antisense ? "GTAG|rev" : "GTAG|fwd", *ref_str, name, cout);
469	>	++itr;
470	>	}
471		}
472
391	–	ref_stream.clear();
392	–	ref_stream.seekg(0, ios::beg);
393	–
394	–
395	–
396	–	while(ref_stream.good() &&
397	–	!ref_stream.eof())
473		{
474	<	Reference ref_str;
475	<	string name;
476	<
477	<	readMeta(ref_stream, name, Fasta());
478	<	string::size_type space_pos = name.find_first_of(" \t\r");
479	<	if (space_pos != string::npos)
480	<	{
481	<	name.resize(space_pos);
482	<	}
408	<
409	<	read(ref_stream, ref_str, Fasta());
410	<
411	<	uint32_t refid = rt.get_id(name, NULL,0);
412	<	Insertion dummy_left(refid, 0, "");
413	<	Insertion dummy_right(refid, VMAXINT32, "");
414	<
415	<	std::set<Insertion>::iterator itr = insertions.lower_bound(dummy_left);
416	<	std::set<Insertion>::iterator upper   = insertions.upper_bound(dummy_right);
474	>	std::set<Deletion>::iterator itr = deletions.begin();
475	>	while(itr != deletions.end())
476	>	{
477	>	RefSequenceTable::Sequence* ref_str = rt.get_seq(itr->refid);
478	>	const char* name = rt.get_name(itr->refid);
479	>	print_splice((Junction)*itr, read_length, itr->antisense ? "del|rev" : "del|fwd", *ref_str, name, cout);
480	>	++itr;
481	>	}
482	>	}
483
484	<	while(itr != upper && itr != insertions.end()){
485	<	print_insertion(*itr, read_length, ref_str, name, cout);
486	<	++itr;
487	<	}
484	>	{
485	>	std::set<Insertion>::iterator itr  = insertions.begin();
486	>	while(itr != insertions.end()){
487	>	RefSequenceTable::Sequence* ref_str = rt.get_seq(itr->refid);
488	>	const char* name = rt.get_name(itr->refid);
489	>	print_insertion(*itr, read_length, *ref_str, name, cout);
490	>	++itr;
491	>	}
492		}
493
494	+	{
495	+	std::set<Fusion>::iterator itr = fusions.begin();
496	+	while(itr != fusions.end()){
497	+	RefSequenceTable::Sequence* left_ref_str = rt.get_seq(itr->refid1);
498	+	RefSequenceTable::Sequence* right_ref_str = rt.get_seq(itr->refid2);
499	+	const char* left_ref_name = rt.get_name(itr->refid1);
500	+	const char* right_ref_name = rt.get_name(itr->refid2);
501	+	print_fusion(*itr, read_length, *left_ref_str, *right_ref_str, left_ref_name, right_ref_name, cout);
502	+	++itr;
503	+	}
504	+	}
505		}
506
507		int main(int argc, char** argv)
#	Line 527 | Line 608
608		}
609
610
611	+	/*
612	+	*/
613	+	string fusion_coords_file_list = argv[optind++];
614	+	vector<string> fusion_coords_file_names;
615	+	vector<FILE*> fusion_coords_files;
616	+	tokenize(fusion_coords_file_list, ",", fusion_coords_file_names);
617	+	for(size_t s = 0; s < fusion_coords_file_names.size(); ++s)
618	+	{
619	+	FILE* fusion_coords_file = fopen(fusion_coords_file_names[s].c_str(),"r");
620	+	if(!fusion_coords_file)
621	+	{
622	+	fprintf(stderr, "Warning: cannot open %s for reading\n",
623	+	fusion_coords_file_names[s].c_str());
624	+	continue;
625	+	}
626	+	fusion_coords_files.push_back(fusion_coords_file);
627	+	}
628	+	if(optind >= argc)
629	+	{
630	+	print_usage();
631	+	return 1;
632	+	}
633	+
634
635		string ref_file_name = argv[optind++];
636		ifstream ref_stream(ref_file_name.c_str());
#	Line 538 | Line 642
642		exit(1);
643		}
644
645	<	driver(coords_files, insertion_coords_files, deletion_coords_files, ref_stream);
645	>	driver(coords_files, insertion_coords_files, deletion_coords_files, fusion_coords_files, ref_stream);
646		return 0;
647		}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines