ViewVC Help

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

root/gclib/tophat_cpp/tophat_reports.cpp

#	Line 30 | Line 30
30		#include <seqan/file.h>
31		#include <getopt.h>
32
33	+	#include <boost/thread.hpp>
34	+	#include <boost/random/mersenne_twister.hpp>
35	+
36		#include "common.h"
37	+	#include "utils.h"
38		#include "bwt_map.h"
39		#include "junctions.h"
40		#include "insertions.h"
41		#include "deletions.h"
42	+	#include "fusions.h"
43		#include "align_status.h"
44		#include "fragments.h"
45		#include "wiggles.h"
46		#include "tokenize.h"
47		#include "reads.h"
48	+	#include "coverage.h"
49	+
50
51		#include "inserts.h"
52
#	Line 47 | Line 54
54		using namespace seqan;
55		using std::set;
56
57	+	static const JunctionSet empty_junctions;
58	+	static const InsertionSet empty_insertions;
59	+	static const DeletionSet empty_deletions;
60	+	static const FusionSet empty_fusions;
61	+	static const Coverage empty_coverage;
62	+
63	+	// daehwan - this is redundancy, which should be removed.
64	+	void get_seqs(istream& ref_stream,
65	+	RefSequenceTable& rt,
66	+	bool keep_seqs = true)
67	+	{
68	+	while(ref_stream.good() && !ref_stream.eof())
69	+	{
70	+	RefSequenceTable::Sequence* ref_str = new RefSequenceTable::Sequence();
71	+	string name;
72	+	readMeta(ref_stream, name, Fasta());
73	+	string::size_type space_pos = name.find_first_of(" \t\r");
74	+	if (space_pos != string::npos)
75	+	{
76	+	name.resize(space_pos);
77	+	}
78	+	fprintf(stderr, "\tLoading %s...", name.c_str());
79	+	seqan::read(ref_stream, *ref_str, Fasta());
80	+	fprintf(stderr, "done\n");
81	+
82	+	rt.get_id(name, keep_seqs ? ref_str : NULL, 0);
83	+	}
84	+	}
85	+
86	+	struct cmp_read_alignment
87	+	{
88	+	bool operator() (const BowtieHit& l, const BowtieHit& r) const
89	+	{
90	+	return l.alignment_score() > r.alignment_score();
91	+	}
92	+	};
93	+
94		void read_best_alignments(const HitsForRead& hits_for_read,
95	<	FragmentAlignmentGrade& best_grade,
96	<	HitsForRead& best_hits,
97	<	const JunctionSet& gtf_junctions)
95	>	HitsForRead& best_hits,
96	>	const JunctionSet& gtf_junctions,
97	>	const JunctionSet& junctions = empty_junctions,
98	>	const InsertionSet& insertions = empty_insertions,
99	>	const DeletionSet& deletions = empty_deletions,
100	>	const FusionSet& fusions = empty_fusions,
101	>	const Coverage& coverage = empty_coverage,
102	>	bool final_report = false)
103		{
104		const vector<BowtieHit>& hits = hits_for_read.hits;
105	+
106	+	if (hits.size() >= max_multihits * 5)
107	+	return;
108	+
109		for (size_t i = 0; i < hits.size(); ++i)
110		{
111	<	if (hits[i].edit_dist()>max_read_mismatches) continue;
112	<	FragmentAlignmentGrade g(hits[i], gtf_junctions);
113	<	// Is the new status better than the current best one?
114	<	if (best_grade < g)
115	<	{
116	<	best_hits.hits.clear();
117	<	best_grade = g;
118	<	best_hits.hits.push_back(hits[i]);
119	<	}
120	<	else if (!(g < best_grade)) // is it just as good?
121	<	{
122	<	best_grade.num_alignments++;
123	<	best_hits.hits.push_back(hits[i]);
124	<	}
111	>	if (hits[i].edit_dist() > max_read_mismatches)
112	>	continue;
113	>
114	>	BowtieHit hit = hits[i];
115	>	AlignStatus align_status(hit, gtf_junctions,
116	>	junctions, insertions, deletions, fusions, coverage);
117	>	hit.alignment_score(align_status._alignment_score);
118	>
119	>	if (report_secondary_alignments || !final_report)
120	>	{
121	>	best_hits.hits.push_back(hit);
122	>	}
123	>	else
124	>	{
125	>	// Is the new status better than the current best one?
126	>	if (best_hits.hits.size() == 0 || cmp_read_alignment()(hit, best_hits.hits[0]))
127	>	{
128	>	best_hits.hits.clear();
129	>	best_hits.hits.push_back(hit);
130	>	}
131	>	else if (!cmp_read_alignment()(best_hits.hits[0], hit)) // is it just as good?
132	>	{
133	>	best_hits.hits.push_back(hit);
134	>	}
135	>	}
136	>	}
137	>
138	>	if ((report_secondary_alignments || !final_report) && best_hits.hits.size() > 0)
139	>	{
140	>	sort(best_hits.hits.begin(), best_hits.hits.end(), cmp_read_alignment());
141	>	}
142	>
143	>	if (final_report)
144	>	{
145	>	if (best_hits.hits.size() > max_multihits)
146	>	best_hits.hits.erase(best_hits.hits.begin() + max_multihits, best_hits.hits.end());
147		}
148		}
149
150	<	void pair_best_alignments(const HitsForRead& left_hits,
76	<	const HitsForRead& right_hits,
77	<	InsertAlignmentGrade& best_grade,
78	<	HitsForRead& left_best_hits,
79	<	HitsForRead& right_best_hits)
150	>	bool set_insert_alignment_grade(const BowtieHit& lh, const BowtieHit& rh, const JunctionSet& junctions, InsertAlignmentGrade& grade)
151		{
152	<	// max mate inner distance (genomic)
153	<	int min_mate_inner_dist = inner_dist_mean - inner_dist_std_dev;
154	<	if (max_mate_inner_dist == -1)
152	>	bool fusion = false;
153	>	bool left_fusion = lh.fusion_opcode() != FUSION_NOTHING;
154	>	bool right_fusion = rh.fusion_opcode() != FUSION_NOTHING;
155	>	if (left_fusion && right_fusion)
156	>	return false;
157	>
158	>	fusion = left_fusion || right_fusion;
159	>	if (!fusion && lh.ref_id() != rh.ref_id())
160	>	fusion = true;
161	>
162	>	if (!fusion && lh.ref_id() == rh.ref_id())
163	>	{
164	>	if (lh.antisense_align() == rh.antisense_align())
165	>	fusion = true;
166	>	else
167		{
168	<	max_mate_inner_dist = inner_dist_mean + inner_dist_std_dev;
168	>	int inner_dist = 0;
169	>	if (lh.antisense_align())
170	>	// used rh.left() instead of rh.right() for the cases,
171	>	// where reads overlap with each other or reads span introns
172	>	inner_dist = lh.left() - rh.left();
173	>	else
174	>	inner_dist = rh.left() - lh.left();
175	>
176	>	if (inner_dist < 0 || inner_dist > (int)fusion_min_dist)
177	>	fusion = true;
178		}
179	+	}
180	+
181	+	// a read contains its partner, in which case the paired mapping will be ignored.
182	+	if (!fusion)
183	+	{
184	+	if (lh.left() <= rh.left() && lh.right() >= rh.right() && lh.right() - lh.left() > rh.right() - rh.left())
185	+	return false;
186	+	else if (rh.left() <= lh.left() && rh.right() >= lh.right() && rh.right() - rh.left() > lh.right() - lh.left())
187	+	return false;
188	+	}
189	+
190	+	grade = InsertAlignmentGrade(lh, rh, junctions, fusion);
191	+
192	+	return true;
193	+	}
194	+
195	+	struct cmp_pair_alignment
196	+	{
197	+	cmp_pair_alignment(const JunctionSet& junctions) :
198	+	_junctions(&junctions) {}
199
200	<	const vector<BowtieHit>& left = left_hits.hits;
201	<	const vector<BowtieHit>& right = right_hits.hits;
202	<
203	<	for (size_t i = 0; i < left.size(); ++i)
200	>	bool operator() (const pair<BowtieHit, BowtieHit>& l, const pair<BowtieHit, BowtieHit>& r) const
201	>	{
202	>	InsertAlignmentGrade gl; set_insert_alignment_grade(l.first, l.second, *_junctions, gl);
203	>	InsertAlignmentGrade gr; set_insert_alignment_grade(r.first, r.second, *_junctions, gr);
204	>
205	>	bool better = gr < gl;
206	>	bool worse = gl < gr;
207	>
208	>	if (better && !worse)
209	>	return true;
210	>	else
211	>	return false;
212	>	}
213	>
214	>	const JunctionSet* _junctions;
215	>	};
216	>
217	>	void pair_best_alignments(const HitsForRead& left_hits,
218	>	const HitsForRead& right_hits,
219	>	InsertAlignmentGrade& best_grade,
220	>	vector<pair<BowtieHit, BowtieHit> >& best_hits,
221	>	const JunctionSet& gtf_junctions,
222	>	const JunctionSet& junctions = empty_junctions,
223	>	const InsertionSet& insertions = empty_insertions,
224	>	const DeletionSet& deletions = empty_deletions,
225	>	const FusionSet& fusions = empty_fusions,
226	>	const Coverage& coverage = empty_coverage,
227	>	bool final_report = false)
228	>	{
229	>	const vector<BowtieHit>& left = left_hits.hits;
230	>	const vector<BowtieHit>& right = right_hits.hits;
231	>
232	>	if (left.size() >= max_multihits * 5 || right.size() >= max_multihits * 5)
233	>	return;
234	>
235	>	for (size_t i = 0; i < left.size(); ++i)
236	>	{
237	>	if (left[i].edit_dist() > max_read_mismatches) continue;
238	>
239	>	BowtieHit lh = left[i];
240	>	AlignStatus align_status(lh, gtf_junctions,
241	>	junctions, insertions, deletions, fusions, coverage);
242	>	lh.alignment_score(align_status._alignment_score);
243	>
244	>	for (size_t j = 0; j < right.size(); ++j)
245		{
246	<	if (left[i].edit_dist()>max_read_mismatches) continue;
246	>	if (right[j].edit_dist() > max_read_mismatches) continue;
247	>
248	>	BowtieHit rh = right[j];
249	>	AlignStatus align_status(rh, gtf_junctions,
250	>	junctions, insertions, deletions, fusions, coverage);
251	>	rh.alignment_score(align_status._alignment_score);
252	>	InsertAlignmentGrade g;
253	>	bool allowed;
254	>	allowed = set_insert_alignment_grade(lh, rh, final_report ? junctions : gtf_junctions, g);
255	>
256	>	// daehwan - for debugging purposes
257	>	#if 0
258	>	if (lh.insert_id() == 325708 && !g.is_fusion() && false)
259	>	{
260	>	fprintf(stderr, "lh %d:%d %s score: %d (from %d) NM: %d\n",
261	>	lh.ref_id(), lh.left(), print_cigar(lh.cigar()).c_str(),
262	>	lh.alignment_score(), left[i].alignment_score(), lh.edit_dist());
263	>	fprintf(stderr, "rh %d:%d %s score: %d (from %d) NM: %d\n",
264	>	rh.ref_id(), rh.left(), print_cigar(rh.cigar()).c_str(),
265	>	rh.alignment_score(), right[j].alignment_score(), rh.edit_dist());
266	>	fprintf(stderr, "combined score: %d is_fusion(%d)\n", g.align_score(), g.is_fusion());
267	>	}
268	>	#endif
269	>
270	>	if (!allowed) continue;
271	>	if (!fusion_search && !report_discordant_pair_alignments && g.is_fusion()) continue;
272
273	<	const BowtieHit& lh = left[i];
274	<	for (size_t j = 0; j < right.size(); ++j)
273	>	if (report_secondary_alignments || !final_report)
274	>	{
275	>	best_hits.push_back(make_pair(lh, rh));
276	>	}
277	>	else
278	>	{
279	>	// Is the new status better than the current best one?
280	>	if (best_grade < g)
281		{
282	<	const BowtieHit& rh = right[j];
283	<
284	<	if (lh.ref_id() != rh.ref_id())
101	<	continue;
102	<	if (right[j].edit_dist()>max_read_mismatches) continue;
103	<	InsertAlignmentGrade g(lh, rh, min_mate_inner_dist, max_mate_inner_dist);
104	<
105	<	// Is the new status better than the current best one?
106	<	if (best_grade < g)
107	<	{
108	<	left_best_hits.hits.clear();
109	<	right_best_hits.hits.clear();
110	<
111	<	best_grade = g;
112	<	left_best_hits.hits.push_back(lh);
113	<	right_best_hits.hits.push_back(rh);
114	<	}
115	<	else if (!(g < best_grade))
116	<	{
117	<	left_best_hits.hits.push_back(lh);
118	<	right_best_hits.hits.push_back(rh);
119	<	best_grade.num_alignments++;
120	<	}
282	>	best_hits.clear();
283	>	best_grade = g;
284	>	best_hits.push_back(make_pair(lh, rh));
285		}
286	+	else if (!(g < best_grade))
287	+	{
288	+	best_hits.push_back(make_pair(lh, rh));
289	+	}
290	+	}
291		}
292	+	}
293	+
294	+	if ((report_secondary_alignments || !final_report) && best_hits.size() > 0)
295	+	{
296	+	cmp_pair_alignment cmp(final_report ? junctions : gtf_junctions);
297	+	sort(best_hits.begin(), best_hits.end(), cmp);
298	+	set_insert_alignment_grade(best_hits[0].first, best_hits[0].second, final_report ? junctions : gtf_junctions, best_grade);
299	+	}
300	+
301	+	if (final_report)
302	+	{
303	+	if (best_hits.size() > max_multihits)
304	+	best_hits.erase(best_hits.begin() + max_multihits, best_hits.end());
305	+	}
306	+
307	+	best_grade.num_alignments = best_hits.size();
308		}
309
310		enum FragmentType {FRAG_UNPAIRED, FRAG_LEFT, FRAG_RIGHT};
#	Line 127 | Line 312
312		void add_auxData(vector<string>& auxdata, vector<string>& sam_toks,
313		const RefSequenceTable& rt,const BowtieHit& bh, FragmentType insert_side,
314		int num_hits, const BowtieHit* next_hit, int hitIndex) {
315	+	bool XS_found = false;
316	+	if (sam_toks.size()>11) {
317
318	<	if (sam_toks.size()>11) {
319	<	for (size_t i=11;i<sam_toks.size();++i) {
320	<	if (sam_toks[i].find("NH:i:")==string::npos)
321	<	auxdata.push_back(sam_toks[i]);
322	<	}
323	<	}
324	<	string aux("NH:i:");
325	<	str_appendInt(aux, num_hits);
318	>	for (size_t i=11;i<sam_toks.size();++i) {
319	>	if (sam_toks[i].find("NH:i:")==string::npos &&
320	>	sam_toks[i].find("XS:i:")==string::npos)
321	>	auxdata.push_back(sam_toks[i]);
322	>
323	>	if (sam_toks[i].find("XS:A:")!=string::npos)
324	>	XS_found = true;
325	>	}
326	>	}
327	>	string aux("NH:i:");
328	>	str_appendInt(aux, num_hits);
329	>	auxdata.push_back(aux);
330	>	if (next_hit) {
331	>	const char* nh_ref_name = "=";
332	>	nh_ref_name = rt.get_name(next_hit->ref_id());
333	>	assert (nh_ref_name != NULL);
334	>	bool same_contig=(next_hit->ref_id()==bh.ref_id());
335	>	aux="CC:Z:";
336	>	aux+= (same_contig ? "=" : nh_ref_name);
337	>	auxdata.push_back(aux);
338	>	aux="CP:i:";
339	>	int nh_gpos=next_hit->left() + 1;
340	>	str_appendInt(aux, nh_gpos);
341		auxdata.push_back(aux);
342	<	if (next_hit) {
141	<	const char* nh_ref_name = "=";
142	<	nh_ref_name = rt.get_name(next_hit->ref_id());
143	<	assert (nh_ref_name != NULL);
144	<	bool same_contig=(next_hit->ref_id()==bh.ref_id());
145	<	aux="CC:Z:";
146	<	aux+= (same_contig ? "=" : nh_ref_name);
147	<	auxdata.push_back(aux);
148	<	aux="CP:i:";
149	<	int nh_gpos=next_hit->left() + 1;
150	<	str_appendInt(aux, nh_gpos);
151	<	auxdata.push_back(aux);
152	<	} //has next_hit
342	>	} //has next_hit
343		// FIXME: this code is still a bit brittle, because it contains no
344		// consistency check that the mates are on opposite strands - a current protocol
345		// requirement, and that the strand indicated by the alignment is consistent
346		// with the orientation of the splices (though that should be handled upstream).
347	+	if (!XS_found) {
348		const string xs_f("XS:A:+");
349		const string xs_r("XS:A:-");
350	<	if (bh.contiguous())  {
351	<	if (library_type == FR_FIRSTSTRAND) {
352	<	if (insert_side == FRAG_LEFT || insert_side == FRAG_UNPAIRED) {
353	<	if (bh.antisense_align())
354	<	auxdata.push_back(xs_f);
355	<	else
356	<	auxdata.push_back(xs_r);
357	<	}
358	<	else {
359	<	if (bh.antisense_align())
360	<	auxdata.push_back(xs_r);
361	<	else
362	<	auxdata.push_back(xs_f);
363	<	}
364	<	}
365	<	else if (library_type == FR_SECONDSTRAND)   {
366	<	if (insert_side == FRAG_LEFT || insert_side == FRAG_UNPAIRED){
367	<	if (bh.antisense_align())
368	<	auxdata.push_back(xs_r);
369	<	else
370	<	auxdata.push_back(xs_f);
371	<	}
372	<	else
373	<	{
374	<	if (bh.antisense_align())
375	<	auxdata.push_back(xs_f);
376	<	else
377	<	auxdata.push_back(xs_r);
378	<	}
379	<	}
380	<	} //bh.contiguous()
190	<	if (hitIndex >= 0)
350	>	if (library_type == FR_FIRSTSTRAND) {
351	>	if (insert_side == FRAG_LEFT || insert_side == FRAG_UNPAIRED) {
352	>	if (bh.antisense_align())
353	>	auxdata.push_back(xs_f);
354	>	else
355	>	auxdata.push_back(xs_r);
356	>	}
357	>	else {
358	>	if (bh.antisense_align())
359	>	auxdata.push_back(xs_r);
360	>	else
361	>	auxdata.push_back(xs_f);
362	>	}
363	>	}
364	>	else if (library_type == FR_SECONDSTRAND)   {
365	>	if (insert_side == FRAG_LEFT || insert_side == FRAG_UNPAIRED){
366	>	if (bh.antisense_align())
367	>	auxdata.push_back(xs_r);
368	>	else
369	>	auxdata.push_back(xs_f);
370	>	}
371	>	else
372	>	{
373	>	if (bh.antisense_align())
374	>	auxdata.push_back(xs_f);
375	>	else
376	>	auxdata.push_back(xs_r);
377	>	}
378	>	}
379	>	}
380	>	if (hitIndex >= 0)
381		{
382	<	string aux("HI:i:");
383	<	str_appendInt(aux, hitIndex);
384	<	auxdata.push_back(aux);
382	>	string aux("HI:i:");
383	>	str_appendInt(aux, hitIndex);
384	>	auxdata.push_back(aux);
385		}
386		}
387
388	<	bool rewrite_sam_record(GBamWriter& bam_writer, const RefSequenceTable& rt,
388	>	bool rewrite_sam_record(GBamWriter& bam_writer,
389	>	const RefSequenceTable& rt,
390		const BowtieHit& bh,
391		const char* bwt_buf,
392		const char* read_alt_name,
202	–	const FragmentAlignmentGrade& grade,
393		FragmentType insert_side,
394		int num_hits,
395		const BowtieHit* next_hit,
396		bool primary,
397		int hitIndex)
398		{
399	<	// Rewrite this hit, filling in the alt name, mate mapping
400	<	// and setting the pair flag
401	<	vector<string> sam_toks;
402	<	tokenize(bwt_buf, "\t", sam_toks);
399	>	// Rewrite this hit, filling in the alt name, mate mapping
400	>	// and setting the pair flag
401	>	vector<string> sam_toks;
402	>	tokenize(bwt_buf, "\t", sam_toks);
403	>
404	>	string ref_name = sam_toks[2], ref_name2 = "";
405	>	char cigar1[1024] = {0}, cigar2[1024] = {0};
406	>	string left_seq, right_seq, left_qual, right_qual;
407	>	int left1 = -1, left2 = -1;
408	>	bool fusion_alignment = false;
409	>	size_t XF_index = 0;
410	>	for (size_t i = 11; i < sam_toks.size(); ++i)
411	>	{
412	>	string& tok = sam_toks[i];
413	>	if (strncmp(tok.c_str(), "XF", 2) == 0)
414	>	{
415	>	fusion_alignment = true;
416	>	XF_index = i;
417	>
418	>	vector<string> tuple_fields;
419	>	tokenize(tok.c_str(), " ", tuple_fields);
420	>	vector<string> contigs;
421	>	tokenize(tuple_fields[1].c_str(), "-", contigs);
422	>	if (contigs.size() >= 2)
423	>	{
424	>	ref_name = contigs[0];
425	>	ref_name2 = contigs[1];
426	>	}
427	>
428	>	extract_partial_hits(bh, tuple_fields[4].c_str(), tuple_fields[5].c_str(),
429	>	cigar1, cigar2, left_seq, right_seq,
430	>	left_qual, right_qual, left1, left2);
431	>
432	>	break;
433	>	}
434	>
435	>	else if (strncmp(tok.c_str(), "AS", 2) == 0)
436	>	{
437	>	char AS_score[128] = {0};
438	>	sprintf(AS_score, "AS:i:%d", min(0, bh.alignment_score()));
439	>	tok = AS_score;
440	>	}
441	>	}
442	>
443		string qname(read_alt_name);
444		size_t slash_pos=qname.rfind('/');
445		if (slash_pos!=string::npos)
446	<	qname.resize(slash_pos);
447	<	//read_alt_name as QNAME
448	<	int flag=atoi(sam_toks[1].c_str()); //FLAG
449	<	if (insert_side != FRAG_UNPAIRED) {
450	<	//flag = atoi(sam_toks[1].c_str());
451	<	// mark this as a singleton mate
452	<	flag |= 0x0001;
453	<	if (insert_side == FRAG_LEFT)
454	<	flag |= 0x0040;
455	<	else if (insert_side == FRAG_RIGHT)
456	<	flag |= 0x0080;
457	<	flag |= 0x0008;
458	<	//char flag_buf[64];
459	<	//sprintf(flag_buf, "%d", flag);
460	<	//sam_toks[t] = flag_buf;
461	<	}
462	<	if (!primary)
463	<	flag |= 0x100;
464	<
465	<	int gpos=isdigit(sam_toks[3][0]) ? atoi(sam_toks[3].c_str()) : 0;
466	<	int mapQ=255;
467	<	if (grade.num_alignments > 1)  {
468	<	double err_prob = 1 - (1.0 / grade.num_alignments);
469	<	mapQ = (int)(-10.0 * log(err_prob) / log(10.0));
470	<	}
471	<	vector<string> auxdata;
472	<	add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex);
473	<	int tlen =atoi(sam_toks[8].c_str()); //TLEN
474	<	int mate_pos=atoi(sam_toks[7].c_str());
475	<	GBamRecord* bamrec=bam_writer.new_record(qname.c_str(), flag, sam_toks[2].c_str(), gpos, mapQ,
476	<	sam_toks[5].c_str(), sam_toks[6].c_str(), mate_pos,
477	<	tlen, sam_toks[9].c_str(), sam_toks[10].c_str(), &auxdata);
478	<	bam_writer.write(bamrec);
479	<	delete bamrec;
480	<	return true;
446	>	qname.resize(slash_pos);
447	>	//read_alt_name as QNAME
448	>	int flag=atoi(sam_toks[1].c_str()); //FLAG
449	>	if (insert_side != FRAG_UNPAIRED) {
450	>	//flag = atoi(sam_toks[1].c_str());
451	>	// mark this as a singleton mate
452	>	flag |= 0x0001;
453	>	if (insert_side == FRAG_LEFT)
454	>	flag |= 0x0040;
455	>	else if (insert_side == FRAG_RIGHT)
456	>	flag |= 0x0080;
457	>	flag |= 0x0008;
458	>	//char flag_buf[64];
459	>	//sprintf(flag_buf, "%d", flag);
460	>	//sam_toks[t] = flag_buf;
461	>	}
462	>	if (!primary)
463	>	flag |= 0x100;
464	>
465	>	int gpos=isdigit(sam_toks[3][0]) ? atoi(sam_toks[3].c_str()) : 0;
466	>	int mapQ = 50;
467	>	if (num_hits > 1)  {
468	>	double err_prob = 1 - (1.0 / num_hits);
469	>	mapQ = (int)(-10.0 * log(err_prob) / log(10.0));
470	>	}
471	>	int tlen =atoi(sam_toks[8].c_str()); //TLEN
472	>	int mate_pos=atoi(sam_toks[7].c_str());
473	>
474	>	string rg_aux = "";
475	>	if (!sam_readgroup_id.empty())
476	>	rg_aux = string("RG:Z:") + sam_readgroup_id;
477	>
478	>	GBamRecord* bamrec=NULL;
479	>	if (fusion_alignment) {
480	>	vector<string> auxdata;
481	>	add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex);
482	>	if (rg_aux != "")
483	>	auxdata.push_back(rg_aux);
484	>	bamrec=bam_writer.new_record(qname.c_str(), flag, ref_name.c_str(), left1 + 1, mapQ,
485	>	cigar1, sam_toks[6].c_str(), mate_pos,
486	>	tlen, left_seq.c_str(), left_qual.c_str(), &auxdata);
487	>	bam_writer.write(bamrec);
488	>	delete bamrec;
489	>
490	>	auxdata.clear();
491	>	sam_toks[XF_index][5] = '2';
492	>	add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex);
493	>	if (rg_aux != "")
494	>	auxdata.push_back(rg_aux);
495	>	bamrec=bam_writer.new_record(qname.c_str(), flag, ref_name2.c_str(), left2 + 1, mapQ,
496	>	cigar2, sam_toks[6].c_str(), mate_pos,
497	>	tlen, right_seq.c_str(), right_qual.c_str(), &auxdata);
498	>	bam_writer.write(bamrec);
499	>	delete bamrec;
500	>	} else {
501	>	vector<string> auxdata;
502	>	add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex);
503	>	if (rg_aux != "")
504	>	auxdata.push_back(rg_aux);
505	>	bamrec=bam_writer.new_record(qname.c_str(), flag, sam_toks[2].c_str(), gpos, mapQ,
506	>	sam_toks[5].c_str(), sam_toks[6].c_str(), mate_pos,
507	>	tlen, sam_toks[9].c_str(), sam_toks[10].c_str(), &auxdata);
508	>	bam_writer.write(bamrec);
509	>	delete bamrec;
510	>	}
511	>
512	>	return true;
513		}
514
515	<	bool rewrite_sam_record(GBamWriter& bam_writer, const RefSequenceTable& rt,
515	>	bool rewrite_sam_record(GBamWriter& bam_writer,
516	>	const RefSequenceTable& rt,
517		const BowtieHit& bh,
518		const char* bwt_buf,
519		const char* read_alt_name,
#	Line 262 | Line 525
525		bool primary,
526		int hitIndex)
527		{
528	<	// Rewrite this hit, filling in the alt name, mate mapping
529	<	// and setting the pair flag
530	<	vector<string> sam_toks;
531	<	tokenize(bwt_buf, "\t", sam_toks);
528	>	// Rewrite this hit, filling in the alt name, mate mapping
529	>	// and setting the pair flag
530	>	vector<string> sam_toks;
531	>	tokenize(bwt_buf, "\t", sam_toks);
532		string qname(read_alt_name);
533		size_t slash_pos=qname.rfind('/');
534		if (slash_pos!=string::npos)
535	<	qname.resize(slash_pos);
536	<	//read_alt_name as QNAME
537	<	int flag = atoi(sam_toks[1].c_str());
538	<	// 0x0010 (strand of query) is assumed to be set correctly
539	<	// to begin with
540	<	flag |= 0x0001; //it must be paired
541	<	if (insert_side == FRAG_LEFT)
542	<	flag |= 0x0040;
543	<	else if (insert_side == FRAG_RIGHT)
544	<	flag |= 0x0080;
545	<	if (!primary)
546	<	flag |= 0x100;
547	<
548	<	int gpos=isdigit(sam_toks[3][0]) ? atoi(sam_toks[3].c_str()) : 0;
549	<	int mapQ=255;
550	<	if (grade.num_alignments > 1) {
551	<	double err_prob = 1 - (1.0 / grade.num_alignments);
552	<	mapQ = (int)(-10.0 * log(err_prob) / log(10.0));
553	<	}
554	<	int tlen=0; //TLEN
555	<	int mate_pos=atoi(sam_toks[7].c_str());
556	<	if (partner) {
557	<	if (partner->ref_id()==bh.ref_id()) {
558	<	sam_toks[6] = "="; //same chromosome
559	<	//TLEN:
560	<	tlen = bh.left() < partner->left() ? partner->right() - bh.left() :
561	<	partner->left() - bh.right();
562	<	}
535	>	qname.resize(slash_pos);
536	>	//read_alt_name as QNAME
537	>	int flag = atoi(sam_toks[1].c_str());
538	>	// 0x0010 (strand of query) is assumed to be set correctly
539	>	// to begin with
540	>	flag |= 0x0001; //it must be paired
541	>	if (insert_side == FRAG_LEFT)
542	>	flag |= 0x0040;
543	>	else if (insert_side == FRAG_RIGHT)
544	>	flag |= 0x0080;
545	>	if (!primary)
546	>	flag |= 0x100;
547	>
548	>	string ref_name = sam_toks[2], ref_name2 = "";
549	>	char cigar1[1024] = {0}, cigar2[1024] = {0};
550	>	string left_seq, right_seq, left_qual, right_qual;
551	>	int left1 = -1, left2 = -1;
552	>	bool fusion_alignment = false;
553	>	size_t XF_tok_idx = 11;
554	>	for (; XF_tok_idx < sam_toks.size(); ++XF_tok_idx)
555	>	{
556	>	string& tok = sam_toks[XF_tok_idx];
557	>	if (strncmp(tok.c_str(), "XF", 2) == 0)
558	>	{
559	>	fusion_alignment = true;
560	>
561	>	vector<string> tuple_fields;
562	>	tokenize(tok.c_str(), " ", tuple_fields);
563	>	vector<string> contigs;
564	>	tokenize(tuple_fields[1].c_str(), "-", contigs);
565	>	if (contigs.size() >= 2)
566	>	{
567	>	ref_name = contigs[0];
568	>	ref_name2 = contigs[1];
569	>	}
570	>
571	>	extract_partial_hits(bh, tuple_fields[4].c_str(), tuple_fields[5].c_str(),
572	>	cigar1, cigar2, left_seq, right_seq,
573	>	left_qual, right_qual, left1, left2);
574	>
575	>	break;
576	>	}
577
578	<	else { //partner on different chromosome/contig
579	<	sam_toks[6] = rt.get_name(partner->ref_id());
580	<	if (sam_toks[6].empty()) {
581	<	//FIXME -- this should never happen
582	<	sam_toks[6] = "=";
583	<	fprintf(stderr, "Warning: partner ref_id %d has no entry in ref table?\n", partner->ref_id());
584	<	}
585	<	}
586	<	mate_pos = partner->left() + 1;
587	<	if (grade.happy())
588	<	flag |= 0x0002;
589	<	if (partner->antisense_align())
590	<	flag |=  0x0020;
591	<	}
592	<	else {
593	<	sam_toks[6] = "*";
594	<	mate_pos = 0;
595	<	flag |= 0x0008;
596	<	}
597	<	vector<string> auxdata;
598	<	add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex);
599	<	GBamRecord* bamrec=bam_writer.new_record(qname.c_str(), flag, sam_toks[2].c_str(), gpos, mapQ,
600	<	sam_toks[5].c_str(), sam_toks[6].c_str(), mate_pos,
601	<	tlen, sam_toks[9].c_str(), sam_toks[10].c_str(), &auxdata);
602	<	bam_writer.write(bamrec);
603	<	delete bamrec;
604	<	return true;
578	>	else if (strncmp(tok.c_str(), "AS", 2) == 0)
579	>	{
580	>	char AS_score[128] = {0};
581	>	sprintf(AS_score, "AS:i:%d", min(0, bh.alignment_score()));
582	>	tok = AS_score;
583	>	}
584	>	}
585	>
586	>	int gpos=isdigit(sam_toks[3][0]) ? atoi(sam_toks[3].c_str()) : 0;
587	>	int mapQ = 50;
588	>	if (grade.num_alignments > 1) {
589	>	double err_prob = 1 - (1.0 / grade.num_alignments);
590	>	mapQ = (int)(-10.0 * log(err_prob) / log(10.0));
591	>	}
592	>	int tlen=0; //TLEN
593	>	int mate_pos=atoi(sam_toks[7].c_str());
594	>	if (partner) {
595	>	if (partner->ref_id()==bh.ref_id()) {
596	>	sam_toks[6] = "="; //same chromosome
597	>	//TLEN:
598	>	tlen = bh.left() < partner->left() ? partner->right() - bh.left() :
599	>	partner->left() - bh.right();
600	>	}
601	>
602	>	else { //partner on different chromosome/contig
603	>	sam_toks[6] = rt.get_name(partner->ref_id());
604	>	}
605	>	mate_pos = partner->left() + 1;
606	>	if (grade.happy())
607	>	flag |= 0x0002;
608	>	if (partner->antisense_align())
609	>	flag |=  0x0020;
610	>
611	>	if (partner->fusion_opcode() != FUSION_NOTHING)
612	>	{
613	>	char partner_pos[1024];
614	>	sprintf(partner_pos, "XP:Z:%s-%s %d", rt.get_name(partner->ref_id()), rt.get_name(partner->ref_id2()), partner->left() + 1);
615	>	sam_toks.push_back(partner_pos);
616	>	}
617	>	}
618	>	else {
619	>	sam_toks[6] = "*";
620	>	mate_pos = 0;
621	>	flag |= 0x0008;
622	>	}
623	>
624	>	string rg_aux = "";
625	>	if (!sam_readgroup_id.empty())
626	>	rg_aux = string("RG:Z:") + sam_readgroup_id;
627	>
628	>	GBamRecord* bamrec=NULL;
629	>	if (fusion_alignment) {
630	>	vector<string> auxdata;
631	>	add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex);
632	>	if (rg_aux != "")
633	>	auxdata.push_back(rg_aux);
634	>	bamrec=bam_writer.new_record(qname.c_str(), flag, ref_name.c_str(), left1 + 1, mapQ,
635	>	cigar1, sam_toks[6].c_str(), mate_pos,
636	>	tlen, left_seq.c_str(), left_qual.c_str(), &auxdata);
637	>	bam_writer.write(bamrec);
638	>	delete bamrec;
639	>
640	>	auxdata.clear();
641	>	sam_toks[XF_tok_idx][5] = '2';
642	>	add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex);
643	>	if (rg_aux != "")
644	>	auxdata.push_back(rg_aux);
645	>	bamrec=bam_writer.new_record(qname.c_str(), flag, ref_name2.c_str(), left2 + 1, mapQ,
646	>	cigar2, sam_toks[6].c_str(), mate_pos,
647	>	tlen, right_seq.c_str(), right_qual.c_str(), &auxdata);
648	>	bam_writer.write(bamrec);
649	>	delete bamrec;
650	>	} else {
651	>	vector<string> auxdata;
652	>	add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex);
653	>	if (rg_aux != "")
654	>	auxdata.push_back(rg_aux);
655	>	bamrec=bam_writer.new_record(qname.c_str(), flag, sam_toks[2].c_str(), gpos, mapQ,
656	>	sam_toks[5].c_str(), sam_toks[6].c_str(), mate_pos,
657	>	tlen, sam_toks[9].c_str(), sam_toks[10].c_str(), &auxdata);
658	>	bam_writer.write(bamrec);
659	>	delete bamrec;
660	>	}
661	>
662	>	return true;
663		}
664
665		struct lex_hit_sort
#	Line 337 | Line 672
672		{
673		const BowtieHit& lhs = _hits.hits[l];
674		const BowtieHit& rhs = _hits.hits[r];
675	<
675	>
676		uint32_t l_id = lhs.ref_id();
677		uint32_t r_id = rhs.ref_id();
678		if (l_id != r_id)
679	<	{
680	<	return (strcmp(_rt.get_name(lhs.ref_id()), _rt.get_name(rhs.ref_id())) < 0);
346	<	}
679	>	return l_id < r_id;
680	>
681		return lhs.left() < rhs.left();
682		}
683
#	Line 352 | Line 686
686		};
687
688		void print_sam_for_single(const RefSequenceTable& rt,
689	<	const HitsForRead& hits,
690	<	const FragmentAlignmentGrade& grade,
691	<	FragmentType frag_type,
692	<	//FILE* reads_file,
693	<	Read& read,
360	<	GBamWriter& bam_writer,
361	<	FILE* um_out//write unmapped reads here
362	<	)
689	>	const HitsForRead& hits,
690	>	FragmentType frag_type,
691	>	const string& alt_name,
692	>	GBamWriter& bam_writer,
693	>	boost::random::mt19937& rng)
694		{
695	<	assert(!read.alt_name.empty());
695	>	//assert(!read.alt_name.empty());
696	>	if (hits.hits.empty())
697	>	return;
698	>
699		lex_hit_sort s(rt, hits);
700		vector<uint32_t> index_vector;
701	<	for (size_t i = 0; i < hits.hits.size(); ++i)
701	>
702	>	size_t i;
703	>	for (i = 0; i < hits.hits.size(); ++i)
704		index_vector.push_back(i);
705
706		sort(index_vector.begin(), index_vector.end(), s);
707	<	size_t primaryHit = (hits.hits.empty()? 0: random() % hits.hits.size());
707	>
708	>	size_t primaryHit = 0;
709	>	if (!report_secondary_alignments)
710	>	primaryHit = rng() % hits.hits.size();
711	>
712		bool multipleHits = (hits.hits.size() > 1);
713	<	for (size_t i = 0; i < hits.hits.size(); ++i)
714	<	{
715	<	bool primary = (i == primaryHit);
716	<	size_t index = index_vector[i];
717	<	const BowtieHit& bh = hits.hits[index];
718	<	rewrite_sam_record(bam_writer, rt,
719	<	bh,
720	<	bh.hitfile_rec().c_str(),
721	<	read.alt_name.c_str(),
722	<	grade,
723	<	frag_type,
724	<	hits.hits.size(),
725	<	(i < hits.hits.size()-1) ? &(hits.hits[index_vector[i+1]]) : NULL,
726	<	primary,
727	<	(multipleHits? i: -1));
388	<	}
713	>	for (i = 0; i < hits.hits.size(); ++i)
714	>	{
715	>	size_t index = index_vector[i];
716	>	bool primary = (index == primaryHit);
717	>	const BowtieHit& bh = hits.hits[index];
718	>	rewrite_sam_record(bam_writer, rt,
719	>	bh,
720	>	bh.hitfile_rec().c_str(),
721	>	alt_name.c_str(),
722	>	frag_type,
723	>	hits.hits.size(),
724	>	(i < hits.hits.size()-1) ? &(hits.hits[index_vector[i+1]]) : NULL,
725	>	primary,
726	>	(multipleHits? i: -1));
727	>	}
728		}
729
730		void print_sam_for_pair(const RefSequenceTable& rt,
731	<	const HitsForRead& left_hits,
393	<	const HitsForRead& right_hits,
731	>	const vector<pair<BowtieHit, BowtieHit> >& best_hits,
732		const InsertAlignmentGrade& grade,
395	–	FILE* left_reads_file,
396	–	FILE* right_reads_file,
733		GBamWriter& bam_writer,
734	<	FILE* left_um_out,
735	<	FILE* right_um_out
736	<	)
734	>	const string& left_alt_name,
735	>	const string& right_alt_name,
736	>	boost::random::mt19937& rng,
737	>	uint64_t begin_id = 0,
738	>	uint64_t end_id = std::numeric_limits<uint64_t>::max())
739		{
402	–	assert (left_hits.insert_id == right_hits.insert_id);
403	–
740		Read left_read;
741		Read right_read;
742	<	assert (left_hits.hits.size() == right_hits.hits.size() ||
743	<	(left_hits.hits.empty() || right_hits.hits.empty()));
742	>	if (best_hits.empty())
743	>	return;
744	>
745	>	size_t i;
746	>	HitsForRead right_hits;
747	>	for (i = 0; i < best_hits.size(); ++i)
748	>	right_hits.hits.push_back(best_hits[i].second);
749
750		size_t primaryHit = 0;
751		vector<uint32_t> index_vector;
752	<	if(right_hits.hits.size() > 0)
753	<	{
754	<	lex_hit_sort s(rt, right_hits);
414	<	for (size_t i = 0; i < right_hits.hits.size(); ++i)
415	<	index_vector.push_back(i);
416	<
417	<	sort(index_vector.begin(), index_vector.end(), s);
418	<	primaryHit = random() % right_hits.hits.size();
419	<	}
420	<	else if (left_hits.hits.size() > 0)
421	<	{
422	<	lex_hit_sort s(rt, left_hits);
423	<	for (size_t i = 0; i < left_hits.hits.size(); ++i)
424	<	index_vector.push_back(i);
425	<
426	<	sort(index_vector.begin(), index_vector.end(), s);
427	<	primaryHit = random() % left_hits.hits.size();
428	<	}
752	>	lex_hit_sort s(rt, right_hits);
753	>	for (i = 0; i < right_hits.hits.size(); ++i)
754	>	index_vector.push_back(i);
755
756	<	bool got_left_read = get_read_from_stream(left_hits.insert_id,
431	<	left_reads_file,
432	<	reads_format,
433	<	false,
434	<	left_read, left_um_out);
435	<
436	<	bool got_right_read = get_read_from_stream(right_hits.insert_id,
437	<	right_reads_file,
438	<	reads_format,
439	<	false,
440	<	right_read, right_um_out);
441	<
442	<	if (left_hits.hits.size() == right_hits.hits.size())
443	<	{
444	<	assert (got_left_read && got_right_read);
445	<	bool multipleHits = (left_hits.hits.size() > 1);
446	<	for (size_t i = 0; i < right_hits.hits.size(); ++i)
447	<	{
448	<	bool primary = (i == primaryHit);
449	<	size_t index = index_vector[i];
450	<	const BowtieHit& right_bh = right_hits.hits[index];
451	<	const BowtieHit& left_bh = left_hits.hits[index];
452	<
453	<	rewrite_sam_record(bam_writer, rt,
454	<	right_bh,
455	<	right_bh.hitfile_rec().c_str(),
456	<	right_read.alt_name.c_str(),
457	<	grade,
458	<	FRAG_RIGHT,
459	<	&left_bh,
460	<	right_hits.hits.size(),
461	<	(i < right_hits.hits.size() - 1) ? &(right_hits.hits[index_vector[i+1]]) : NULL,
462	<	primary,
463	<	(multipleHits? i: -1));
464	<	rewrite_sam_record(bam_writer, rt,
465	<	left_bh,
466	<	left_bh.hitfile_rec().c_str(),
467	<	left_read.alt_name.c_str(),
468	<	grade,
469	<	FRAG_LEFT,
470	<	&right_bh,
471	<	left_hits.hits.size(),
472	<	(i < left_hits.hits.size() - 1) ? &(left_hits.hits[index_vector[i+1]]) : NULL,
473	<	primary,
474	<	(multipleHits? i: -1));
475	<	}
476	<	}
477	<	else if (left_hits.hits.empty())
478	<	{ //only right read was mapped properly
479	<	if (right_um_out) {
480	<	//write it in the mapped file with the #MAPPED# flag
481	<	fprintf(right_um_out, "@%s #MAPPED#\n%s\n+\n%s\n", right_read.alt_name.c_str(),
482	<	right_read.seq.c_str(), right_read.qual.c_str());
483	<	}
484	<	for (size_t i = 0; i < right_hits.hits.size(); ++i)
485	<	{
486	<	bool primary = (i == primaryHit);
487	<	size_t index = index_vector[i];
488	<	const BowtieHit& bh = right_hits.hits[index];
489	<
490	<	rewrite_sam_record(bam_writer, rt,
491	<	bh,
492	<	bh.hitfile_rec().c_str(),
493	<	right_read.alt_name.c_str(),
494	<	grade,
495	<	FRAG_RIGHT,
496	<	NULL,
497	<	right_hits.hits.size(),
498	<	(i < right_hits.hits.size() - 1) ? &(right_hits.hits[index_vector[i+1]]) : NULL,
499	<	primary,
500	<	-1);
501	<
502	<	}
503	<	}
504	<	else if (right_hits.hits.empty())
505	<	{ //only left read was mapped properly
506	<	if (left_um_out) {
507	<	//write it in the mapped file with the #MAPPED# flag
508	<	fprintf(left_um_out, "@%s #MAPPED#\n%s\n+\n%s\n", left_read.alt_name.c_str(), left_read.seq.c_str(),
509	<	left_read.qual.c_str());
510	<	}
756	>	sort(index_vector.begin(), index_vector.end(), s);
757
758	<	for (size_t i = 0; i < left_hits.hits.size(); ++i)
759	<	{
760	<	bool primary = (i == primaryHit);
761	<	size_t index = index_vector[i];
762	<	const BowtieHit& bh = left_hits.hits[index];
763	<	rewrite_sam_record(bam_writer, rt,
764	<	bh,
765	<	bh.hitfile_rec().c_str(),
766	<	left_read.alt_name.c_str(),
767	<	grade,
768	<	FRAG_LEFT,
769	<	NULL,
770	<	left_hits.hits.size(),
771	<	(i < left_hits.hits.size() - 1) ? &(left_hits.hits[index_vector[i+1]]) : NULL,
772	<	primary,
773	<	-1);
774	<
775	<	}
776	<	}
777	<	else
778	<	{
779	<	assert (false);
780	<	}
758	>	if (!report_secondary_alignments)
759	>	primaryHit = rng() % right_hits.hits.size();
760	>
761	>	bool multipleHits = (best_hits.size() > 1);
762	>	for (i = 0; i < best_hits.size(); ++i)
763	>	{
764	>	size_t index = index_vector[i];
765	>	bool primary = (index == primaryHit);
766	>	const BowtieHit& right_bh = best_hits[index].second;
767	>	const BowtieHit& left_bh = best_hits[index].first;
768	>
769	>	rewrite_sam_record(bam_writer, rt,
770	>	right_bh,
771	>	right_bh.hitfile_rec().c_str(),
772	>	right_alt_name.c_str(),
773	>	grade,
774	>	FRAG_RIGHT,
775	>	&left_bh,
776	>	best_hits.size(),
777	>	(i < best_hits.size() - 1) ? &(best_hits[index_vector[i+1]].second) : NULL,
778	>	primary,
779	>	(multipleHits? i: -1));
780	>	rewrite_sam_record(bam_writer, rt,
781	>	left_bh,
782	>	left_bh.hitfile_rec().c_str(),
783	>	left_alt_name.c_str(),
784	>	grade,
785	>	FRAG_LEFT,
786	>	&right_bh,
787	>	best_hits.size(),
788	>	(i < best_hits.size() - 1) ? &(best_hits[index_vector[i+1]].first) : NULL,
789	>	primary,
790	>	(multipleHits? i: -1));
791	>	}
792		}
793
794		/**
#	Line 544 | Line 801
801		InsertionSet& insertions,
802		DeletionSet& deletions)
803		{
804	<	for (size_t i = 0; i < hits.hits.size(); ++i)
804	>	for (size_t i = 0; i < hits.hits.size(); ++i)
805	>	{
806	>	const BowtieHit& bh = hits.hits[i];
807	>	insertions_from_alignment(bh, insertions);
808	>	deletions_from_alignment(bh, deletions);
809	>	}
810	>	}
811	>
812	>	void update_coverage(const HitsForRead& hits,
813	>	Coverage& coverage)
814	>	{
815	>	for (size_t i = 0; i < hits.hits.size(); ++i)
816	>	{
817	>	const BowtieHit& hit = hits.hits[i];
818	>	const vector<CigarOp>& cigar = hit.cigar();
819	>	unsigned int positionInGenome = hit.left();
820	>	RefID ref_id = hit.ref_id();
821	>
822	>	for(size_t c = 0; c < cigar.size(); ++c)
823		{
824	<	const BowtieHit& bh = hits.hits[i];
825	<	insertions_from_alignment(bh, insertions);
826	<	deletions_from_alignment(bh, deletions);
824	>	int opcode = cigar[c].opcode;
825	>	int length = cigar[c].length;
826	>	switch(opcode)
827	>	{
828	>	case REF_SKIP:
829	>	case MATCH:
830	>	case DEL:
831	>	if (opcode == MATCH)
832	>	coverage.add_coverage(ref_id, positionInGenome, length);
833	>
834	>	positionInGenome += length;
835	>	break;
836	>	case rEF_SKIP:
837	>	case mATCH:
838	>	case dEL:
839	>	positionInGenome -= length;
840	>
841	>	if (opcode == mATCH)
842	>	coverage.add_coverage(ref_id, positionInGenome + 1, length);
843	>	break;
844	>	case FUSION_FF:
845	>	case FUSION_FR:
846	>	case FUSION_RF:
847	>	case FUSION_RR:
848	>	positionInGenome = length;
849	>	ref_id = hit.ref_id2();
850	>	break;
851	>	default:
852	>	break;
853	>	}
854		}
855	+	}
856		}
857
858
859	+	void update_fusions(const HitsForRead& hits,
860	+	RefSequenceTable& rt,
861	+	FusionSet& fusions,
862	+	const FusionSet& fusions_ref = empty_fusions)
863	+	{
864	+	if (!fusion_search)
865	+	return;
866	+
867	+	if (hits.hits.size() > fusion_multireads)
868	+	return;
869	+
870	+	bool update_stat = fusions_ref.size() > 0;
871	+	for (size_t i = 0; i < hits.hits.size(); ++i)
872	+	{
873	+	const BowtieHit& bh = hits.hits[i];
874	+
875	+	if (bh.edit_dist() > fusion_read_mismatches)
876	+	continue;
877	+
878	+	fusions_from_alignment(bh, fusions, rt, update_stat);
879	+
880	+	if (update_stat)
881	+	unsupport_fusions(bh, fusions, fusions_ref);
882	+	}
883	+	}
884
885		void update_junctions(const HitsForRead& hits,
886	<	JunctionSet& junctions)
886	>	JunctionSet& junctions)
887		{
888	<	for (size_t i = 0; i < hits.hits.size(); ++i)
889	<	{
890	<	const BowtieHit& bh = hits.hits[i];
891	<	junctions_from_alignment(bh, junctions);
892	<	}
888	>	for (size_t i = 0; i < hits.hits.size(); ++i)
889	>	{
890	>	const BowtieHit& bh = hits.hits[i];
891	>	junctions_from_alignment(bh, junctions);
892	>	}
893		}
894
895		// Extracts junctions from all the SAM hits (based on REF_SKIPs) in the hit file
896		// resets the stream when finished.
897	<	void exclude_hits_on_filtered_junctions(const JunctionSet& junctions,
570	<	HitsForRead& hits)
897	>	void exclude_hits_on_filtered_junctions(const JunctionSet& junctions, HitsForRead& hits)
898		{
899	<	HitsForRead remaining;
900	<	remaining.insert_id = hits.insert_id;
901	<
902	<	for (size_t i = 0; i < hits.hits.size(); ++i)
899	>	HitsForRead remaining;
900	>	remaining.insert_id = hits.insert_id;
901	>
902	>	for (size_t i = 0; i < hits.hits.size(); ++i)
903	>	{
904	>	BowtieHit& bh = hits.hits[i];
905	>	if (bh.edit_dist() > max_read_mismatches)
906	>	continue;
907	>
908	>	bool filter_hit = false;
909	>	if (!bh.contiguous())
910		{
911	<	BowtieHit& bh = hits.hits[i];
912	<	bool filter_hit = false;
913	<	if (!bh.contiguous())
914	<	{
915	<	JunctionSet bh_juncs;
916	<	junctions_from_alignment(bh, bh_juncs);
917	<	for (JunctionSet::iterator itr = bh_juncs.begin();
918	<	itr != bh_juncs.end();
919	<	itr++)
586	<	{
587	<	const Junction& j = itr->first;
588	<	JunctionSet::const_iterator target = junctions.find(j);
589	<	if (target == junctions.end() || !target->second.accepted)
590	<	{
591	<	filter_hit = true;
592	<	break;
593	<	}
594	<	}
595	<	}
596	<	if (!filter_hit)
911	>	JunctionSet bh_juncs;
912	>	junctions_from_alignment(bh, bh_juncs);
913	>	for (JunctionSet::iterator itr = bh_juncs.begin();
914	>	itr != bh_juncs.end();
915	>	itr++)
916	>	{
917	>	const Junction& j = itr->first;
918	>	JunctionSet::const_iterator target = junctions.find(j);
919	>	if (target == junctions.end() || !target->second.accepted)
920		{
921	<	remaining.hits.push_back(bh);
921	>	filter_hit = true;
922	>	break;
923		}
924	<	else
924	>	}
925	>	}
926	>	if (!filter_hit)
927	>	remaining.hits.push_back(bh);
928	>	}
929	>	hits = remaining;
930	>	}
931	>
932	>	void realign_reads(HitsForRead& hits,
933	>	const RefSequenceTable& rt,
934	>	const JunctionSet& junctions,
935	>	const JunctionSet& rev_junctions,
936	>	const InsertionSet& insertions,
937	>	const DeletionSet& deletions,
938	>	const DeletionSet& rev_deletions,
939	>	const FusionSet& fusions)
940	>	{
941	>	if (color)
942	>	return;
943	>
944	>	vector<BowtieHit> additional_hits;
945	>	for (size_t i = 0; i < hits.hits.size(); ++i)
946	>	{
947	>	BowtieHit& bh = hits.hits[i];
948	>	if (fusion_search && bh.fusion_opcode() != FUSION_NOTHING)
949	>	return;
950	>
951	>	const vector<CigarOp>& cigars = bh.cigar();
952	>	int pos = bh.left();
953	>	int refid = bh.ref_id();
954	>	for (size_t j = 0; j < cigars.size(); ++j)
955	>	{
956	>	const CigarOp& op = cigars[j];
957	>	if (j == 0 || j == cigars.size() - 1)
958	>	{
959	>	// let's do this for MATCH case only,
960	>	if (op.opcode == MATCH || op.opcode == mATCH)
961		{
962	<	//fprintf(stderr, "Filtering hit\n");
962	>	int left1, left2;
963	>	if (op.opcode == MATCH)
964	>	{
965	>	left1 = pos;
966	>	left2 = pos + op.length - 1;
967	>	}
968	>	else
969	>	{
970	>	left1 = pos - op.length + 1;
971	>	left2 = pos;
972	>	}
973	>
974	>	{
975	>	JunctionSet::const_iterator lb, ub;
976	>	JunctionSet temp_junctions;
977	>	if (j == 0)
978	>	{
979	>	lb = rev_junctions.upper_bound(Junction(refid, left1, 0, true));
980	>	ub = rev_junctions.lower_bound(Junction(refid, left2, left2, false));
981	>	while (lb != ub && lb != rev_junctions.end())
982	>	{
983	>	Junction temp_junction = lb->first;
984	>	temp_junction.left = lb->first.right;
985	>	temp_junction.right = lb->first.left;
986	>	temp_junctions[temp_junction] = lb->second;
987	>	++lb;
988	>	}
989	>	}
990	>
991	>	if (j == cigars.size() - 1)
992	>	{
993	>	int common_right = left2 + max_report_intron_length;
994	>	lb = junctions.upper_bound(Junction(refid, left1, common_right, true));
995	>	ub = junctions.lower_bound(Junction(refid, left2, common_right, false));
996	>	while (lb != ub && lb != junctions.end())
997	>	{
998	>	temp_junctions[lb->first] = lb->second;
999	>	++lb;
1000	>	}
1001	>	}
1002	>
1003	>	if (temp_junctions.size() > 10)
1004	>	continue;
1005	>
1006	>	JunctionSet::const_iterator junc_iter = temp_junctions.begin();
1007	>	for (; junc_iter != temp_junctions.end(); ++junc_iter)
1008	>	{
1009	>	Junction junc = junc_iter->first;
1010	>
1011	>	/*
1012	>	fprintf(stderr, "%d %d-%d %s (AS:%d XM:%d) with junc %u-%u\n",
1013	>	bh.insert_id(), bh.left(), bh.right(),
1014	>	print_cigar(bh.cigar()).c_str(),
1015	>	bh.alignment_score(), bh.edit_dist(),
1016	>	junc.left, junc.right);
1017	>	//*/
1018	>
1019	>	int new_left = bh.left();
1020	>	int intron_length = junc.right - junc.left - 1;
1021	>	vector<CigarOp> new_cigars;
1022	>	bool anchored = false;
1023	>	if (j == 0 && bh.left() > (int)junc.left)
1024	>	{
1025	>	new_left -= intron_length;
1026	>	int before_match_length = junc.left - new_left + 1;;
1027	>	int after_match_length = op.length - before_match_length;
1028	>
1029	>	if (before_match_length > 0 && after_match_length)
1030	>	{
1031	>	anchored = true;
1032	>	new_cigars.push_back(CigarOp(MATCH, before_match_length));
1033	>	new_cigars.push_back(CigarOp(REF_SKIP, intron_length));
1034	>	new_cigars.push_back(CigarOp(MATCH, after_match_length));
1035	>
1036	>	new_cigars.insert(new_cigars.end(), cigars.begin() + 1, cigars.end());
1037	>	}
1038	>	}
1039	>	else
1040	>	{
1041	>	new_cigars.insert(new_cigars.end(), cigars.begin(), cigars.end() - 1);
1042	>
1043	>	int before_match_length = junc.left - pos + 1;
1044	>	int after_match_length = op.length - before_match_length;
1045	>
1046	>	if (before_match_length > 0 && after_match_length > 0)
1047	>	{
1048	>	anchored = true;
1049	>	new_cigars.push_back(CigarOp(MATCH, before_match_length));
1050	>	new_cigars.push_back(CigarOp(REF_SKIP, intron_length));
1051	>	new_cigars.push_back(CigarOp(MATCH, after_match_length));
1052	>	}
1053	>	}
1054	>
1055	>	BowtieHit new_bh(bh.ref_id(),
1056	>	bh.ref_id2(),
1057	>	bh.insert_id(),
1058	>	new_left,
1059	>	new_cigars,
1060	>	bh.antisense_align(),
1061	>	junc.antisense,
1062	>	0, /* edit_dist - needs to be recalculated */
1063	>	0, /* splice_mms - needs to be recalculated */
1064	>	false);
1065	>
1066	>	new_bh.seq(bh.seq());
1067	>	new_bh.qual(bh.qual());
1068	>
1069	>	const RefSequenceTable::Sequence* ref_str = rt.get_seq(bh.ref_id());
1070	>
1071	>	if (new_left >= 0 &&
1072	>	new_bh.right() <= (int)length(*ref_str) &&
1073	>	anchored)
1074	>	{
1075	>	vector<string> aux_fields;
1076	>	bowtie_sam_extra(new_bh, rt, aux_fields);
1077	>
1078	>	vector<string>::const_iterator aux_iter = aux_fields.begin();
1079	>	for (; aux_iter != aux_fields.end(); ++aux_iter)
1080	>	{
1081	>	const string& aux_field = *aux_iter;
1082	>	if (strncmp(aux_field.c_str(), "AS", 2) == 0)
1083	>	{
1084	>	int alignment_score = atoi(aux_field.c_str() + 5);
1085	>	new_bh.alignment_score(alignment_score);
1086	>	}
1087	>	else if (strncmp(aux_field.c_str(), "XM", 2) == 0)
1088	>	{
1089	>	int XM_value =  atoi(aux_field.c_str() + 5);
1090	>	new_bh.edit_dist(XM_value);
1091	>	}
1092	>	}
1093	>
1094	>	vector<string> sam_toks;
1095	>	tokenize(bh.hitfile_rec().c_str(), "\t", sam_toks);
1096	>
1097	>	char coord[20] = {0,};
1098	>	sprintf(coord, "%d", new_bh.left() + 1);
1099	>	sam_toks[3] = coord;
1100	>	sam_toks[5] = print_cigar(new_bh.cigar());
1101	>	for (size_t a = 11; a < sam_toks.size(); ++a)
1102	>	{
1103	>	string& sam_tok = sam_toks[a];
1104	>	for (size_t b = 0; b < aux_fields.size(); ++b)
1105	>	{
1106	>	const string& aux_tok = aux_fields[b];
1107	>	if (strncmp(sam_tok.c_str(), aux_tok.c_str(), 5) == 0)
1108	>	{
1109	>	sam_tok = aux_tok;
1110	>	break;
1111	>	}
1112	>	}
1113	>	}
1114	>
1115	>	if (!bh.is_spliced())
1116	>	{
1117	>	if (junc.antisense)
1118	>	sam_toks.push_back("XS:A:-");
1119	>	else
1120	>	sam_toks.push_back("XS:A:+");
1121	>	}
1122	>
1123	>
1124	>	string new_rec = "";
1125	>	for (size_t d = 0; d < sam_toks.size(); ++d)
1126	>	{
1127	>	new_rec += sam_toks[d];
1128	>	if (d < sam_toks.size() - 1)
1129	>	new_rec += "\t";
1130	>	}
1131	>
1132	>	new_bh.hitfile_rec(new_rec);
1133	>
1134	>	if (new_bh.edit_dist() <= bh.edit_dist())
1135	>	additional_hits.push_back(new_bh);
1136	>
1137	>	/*
1138	>	fprintf(stderr, "\t%d %d-%d %s (AS:%d XM:%d) with junc %u-%u\n",
1139	>	new_bh.insert_id(), new_bh.left(), new_bh.right(),
1140	>	print_cigar(new_bh.cigar()).c_str(),
1141	>	new_bh.alignment_score(), new_bh.edit_dist(),
1142	>	junc.left, junc.right);
1143	>	//*/
1144	>	}
1145	>	}
1146	>	}
1147	>
1148	>	#if 0
1149	>	{
1150	>	DeletionSet::const_iterator lb, ub;
1151	>	bool use_rev_deletions = (j == 0);
1152	>	const DeletionSet& curr_deletions = (use_rev_deletions ? rev_deletions : deletions);
1153	>	if (use_rev_deletions)
1154	>	{
1155	>	lb = curr_deletions.upper_bound(Deletion(refid, left1, 0, true));
1156	>	ub = curr_deletions.lower_bound(Deletion(refid, left2, left2, false));
1157	>	}
1158	>	else
1159	>	{
1160	>	int common_right = left2 + 100;
1161	>	lb = curr_deletions.upper_bound(Deletion(refid, left1, common_right, true));
1162	>	ub = curr_deletions.lower_bound(Deletion(refid, left2, common_right, false));
1163	>	}
1164	>
1165	>	while (lb != curr_deletions.end() && lb != ub)
1166	>	{
1167	>	Deletion del = lb->first;
1168	>	if (use_rev_deletions)
1169	>	{
1170	>	int temp = del.left;
1171	>	del.left = del.right;
1172	>	del.right = temp;
1173	>	}
1174	>
1175	>	// daehwan - for debuggin purposes
1176	>	/*
1177	>	fprintf(stderr, "(type%d) %d %d-%d %s (AS:%d XM:%d) with junc %u-%u\n",
1178	>	!use_rev_junctions,
1179	>	bh.insert_id(), bh.left(), bh.right(),
1180	>	print_cigar(bh.cigar()).c_str(),
1181	>	bh.alignment_score(), bh.edit_dist(),
1182	>	junc.left, junc.right);
1183	>	*/
1184	>
1185	>	int del_length = del.right - del.left - 1;
1186	>	int new_left = bh.left();
1187	>	if (j == 0)
1188	>	new_left -= del_length;
1189	>
1190	>	vector<CigarOp> new_cigars;
1191	>	if (j == 0)
1192	>	{
1193	>	int before_match_length = del.left - new_left + 1;;
1194	>	int after_match_length = op.length - before_match_length;
1195	>
1196	>	if (before_match_length > 0)
1197	>	new_cigars.push_back(CigarOp(MATCH, before_match_length));
1198	>	new_cigars.push_back(CigarOp(DEL, del_length));
1199	>	if (after_match_length > 0)
1200	>	new_cigars.push_back(CigarOp(MATCH, after_match_length));
1201	>
1202	>	new_cigars.insert(new_cigars.end(), cigars.begin() + 1, cigars.end());
1203	>	}
1204	>	else
1205	>	{
1206	>	new_cigars.insert(new_cigars.end(), cigars.begin(), cigars.end() - 1);
1207	>
1208	>	int before_match_length = del.left - pos + 1;
1209	>	int after_match_length = op.length - before_match_length;
1210	>
1211	>	if (before_match_length > 0)
1212	>	new_cigars.push_back(CigarOp(MATCH, before_match_length));
1213	>	new_cigars.push_back(CigarOp(DEL, del_length));
1214	>	if (after_match_length > 0)
1215	>	new_cigars.push_back(CigarOp(MATCH, after_match_length));
1216	>	}
1217	>
1218	>	BowtieHit new_bh(bh.ref_id(),
1219	>	bh.ref_id2(),
1220	>	bh.insert_id(),
1221	>	new_left,
1222	>	new_cigars,
1223	>	bh.antisense_align(),
1224	>	bh.antisense_splice(),
1225	>	0, /* edit_dist - needs to be recalculated */
1226	>	0, /* splice_mms - needs to be recalculated */
1227	>	false);
1228	>
1229	>	new_bh.seq(bh.seq());
1230	>	new_bh.qual(bh.qual());
1231	>
1232	>	vector<string> aux_fields;
1233	>	bowtie_sam_extra(new_bh, rt, aux_fields);
1234	>
1235	>	vector<string>::const_iterator aux_iter = aux_fields.begin();
1236	>	for (; aux_iter != aux_fields.end(); ++aux_iter)
1237	>	{
1238	>	const string& aux_field = *aux_iter;
1239	>	if (strncmp(aux_field.c_str(), "AS", 2) == 0)
1240	>	{
1241	>	int alignment_score = atoi(aux_field.c_str() + 5);
1242	>	new_bh.alignment_score(alignment_score);
1243	>	}
1244	>	else if (strncmp(aux_field.c_str(), "XM", 2) == 0)
1245	>	{
1246	>	int XM_value =  atoi(aux_field.c_str() + 5);
1247	>	new_bh.edit_dist(XM_value);
1248	>	}
1249	>	}
1250	>
1251	>	vector<string> sam_toks;
1252	>	tokenize(bh.hitfile_rec().c_str(), "\t", sam_toks);
1253	>
1254	>	char coord[20] = {0,};
1255	>	sprintf(coord, "%d", new_bh.left() + 1);
1256	>	sam_toks[3] = coord;
1257	>	sam_toks[5] = print_cigar(new_bh.cigar());
1258	>	for (size_t a = 11; a < sam_toks.size(); ++a)
1259	>	{
1260	>	string& sam_tok = sam_toks[a];
1261	>	for (size_t b = 0; b < aux_fields.size(); ++b)
1262	>	{
1263	>	const string& aux_tok = aux_fields[b];
1264	>	if (strncmp(sam_tok.c_str(), aux_tok.c_str(), 5) == 0)
1265	>	{
1266	>	sam_tok = aux_tok;
1267	>	break;
1268	>	}
1269	>	}
1270	>	}
1271	>
1272	>	string new_rec = "";
1273	>	for (size_t d = 0; d < sam_toks.size(); ++d)
1274	>	{
1275	>	new_rec += sam_toks[d];
1276	>	if (d < sam_toks.size() - 1)
1277	>	new_rec += "\t";
1278	>	}
1279	>
1280	>	new_bh.hitfile_rec(new_rec);
1281	>
1282	>	if (new_bh.edit_dist() <= bh.edit_dist())
1283	>	additional_hits.push_back(new_bh);
1284	>
1285	>	/*
1286	>	fprintf(stderr, "\t%d %d-%d %s (AS:%d XM:%d) with junc %u-%u\n",
1287	>	new_bh.insert_id(), new_bh.left(), new_bh.right(),
1288	>	print_cigar(new_bh.cigar()).c_str(),
1289	>	new_bh.alignment_score(), new_bh.edit_dist(),
1290	>	junc.left, junc.right);
1291	>	*/
1292	>
1293	>	++lb;
1294	>	}
1295	>	}
1296	>
1297	>	{
1298	>	InsertionSet::const_iterator lb, ub;
1299	>	lb = insertions.upper_bound(Insertion(refid, left1, ""));
1300	>	ub = insertions.lower_bound(Insertion(refid, left2, ""));
1301	>
1302	>	while (lb != insertions.end() && lb != ub)
1303	>	{
1304	>	// daehwan - for debugging purposse
1305	>	break;
1306	>
1307	>	Insertion ins = lb->first;
1308	>
1309	>	// daehwan - for debugging purposes
1310	>	/*
1311	>	fprintf(stderr, "(type%d) %d %d-%d %s (AS:%d XM:%d) with junc %u-%u\n",
1312	>	!use_rev_junctions,
1313	>	bh.insert_id(), bh.left(), bh.right(),
1314	>	print_cigar(bh.cigar()).c_str(),
1315	>	bh.alignment_score(), bh.edit_dist(),
1316	>	junc.left, junc.right);
1317	>	*/
1318	>
1319	>	int ins_length = ins.sequence.length();
1320	>	int new_left = bh.left();
1321	>	if (j == 0)
1322	>	new_left -= ins_length;
1323	>
1324	>	vector<CigarOp> new_cigars;
1325	>	if (j == 0)
1326	>	{
1327	>	int before_match_length = ins.left - new_left + 1;;
1328	>	int after_match_length = op.length - before_match_length - ins_length;
1329	>
1330	>	if (before_match_length > 0)
1331	>	new_cigars.push_back(CigarOp(MATCH, before_match_length));
1332	>	new_cigars.push_back(CigarOp(INS, ins_length));
1333	>	if (after_match_length > 0)
1334	>	new_cigars.push_back(CigarOp(MATCH, after_match_length));
1335	>
1336	>	new_cigars.insert(new_cigars.end(), cigars.begin() + 1, cigars.end());
1337	>	}
1338	>	else
1339	>	{
1340	>	new_cigars.insert(new_cigars.end(), cigars.begin(), cigars.end() - 1);
1341	>
1342	>	int before_match_length = ins.left - pos + 1;
1343	>	int after_match_length = op.length - before_match_length - ins_length;
1344	>
1345	>	if (before_match_length > 0)
1346	>	new_cigars.push_back(CigarOp(MATCH, before_match_length));
1347	>	new_cigars.push_back(CigarOp(INS, ins_length));
1348	>	if (after_match_length > 0)
1349	>	new_cigars.push_back(CigarOp(MATCH, after_match_length));
1350	>	}
1351	>
1352	>	BowtieHit new_bh(bh.ref_id(),
1353	>	bh.ref_id2(),
1354	>	bh.insert_id(),
1355	>	new_left,
1356	>	new_cigars,
1357	>	bh.antisense_align(),
1358	>	bh.antisense_splice(),
1359	>	0, /* edit_dist - needs to be recalculated */
1360	>	0, /* splice_mms - needs to be recalculated */
1361	>	false);
1362	>
1363	>	new_bh.seq(bh.seq());
1364	>	new_bh.qual(bh.qual());
1365	>
1366	>	vector<string> aux_fields;
1367	>	bowtie_sam_extra(new_bh, rt, aux_fields);
1368	>
1369	>	vector<string>::const_iterator aux_iter = aux_fields.begin();
1370	>	for (; aux_iter != aux_fields.end(); ++aux_iter)
1371	>	{
1372	>	const string& aux_field = *aux_iter;
1373	>	if (strncmp(aux_field.c_str(), "AS", 2) == 0)
1374	>	{
1375	>	int alignment_score = atoi(aux_field.c_str() + 5);
1376	>	new_bh.alignment_score(alignment_score);
1377	>	}
1378	>	else if (strncmp(aux_field.c_str(), "XM", 2) == 0)
1379	>	{
1380	>	int XM_value =  atoi(aux_field.c_str() + 5);
1381	>	new_bh.edit_dist(XM_value);
1382	>	}
1383	>	}
1384	>
1385	>	/*
1386	>	fprintf(stderr, "\t%d %d-%d %s (AS:%d XM:%d) with junc %u-%u\n",
1387	>	new_bh.insert_id(), new_bh.left(), new_bh.right(),
1388	>	print_cigar(new_bh.cigar()).c_str(),
1389	>	new_bh.alignment_score(), new_bh.edit_dist(),
1390	>	junc.left, junc.right);
1391	>	*/
1392	>
1393	>	++lb;
1394	>	}
1395	>	}
1396	>	#endif
1397		}
1398	+	}
1399	+
1400	+	switch(op.opcode)
1401	+	{
1402	+	case REF_SKIP:
1403	+	pos += op.length;
1404	+	break;
1405	+	case rEF_SKIP:
1406	+	pos -= op.length;
1407	+	break;
1408	+	case MATCH:
1409	+	case DEL:
1410	+	pos += op.length;
1411	+	break;
1412	+	case mATCH:
1413	+	case dEL:
1414	+	pos -= op.length;
1415	+	break;
1416	+	case FUSION_FF:
1417	+	case FUSION_FR:
1418	+	case FUSION_RF:
1419	+	case FUSION_RR:
1420	+	pos = op.length;
1421	+	refid = bh.ref_id2();
1422	+	break;
1423	+	default:
1424	+	break;
1425	+	}
1426		}
1427	<	hits = remaining;
1427	>	}
1428	>
1429	>	hits.hits.insert(hits.hits.end(), additional_hits.begin(), additional_hits.end());
1430	>
1431	>	std::sort(hits.hits.begin(), hits.hits.end());
1432	>	vector<BowtieHit>::iterator new_end = std::unique(hits.hits.begin(), hits.hits.end());
1433	>	hits.hits.erase(new_end, hits.hits.end());
1434		}
1435
1436	<	void get_junctions_from_best_hits(HitStream& left_hs,
1437	<	HitStream& right_hs,
1438	<	ReadTable& it,
1439	<	JunctionSet& junctions,
1440	<	const JunctionSet& gtf_junctions)
1441	<	{
1442	<	HitsForRead curr_left_hit_group;
1443	<	HitsForRead curr_right_hit_group;
1444	<
1445	<	left_hs.next_read_hits(curr_left_hit_group);
1446	<	right_hs.next_read_hits(curr_right_hit_group);
1447	<
1448	<	uint32_t curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id);
1449	<	uint32_t curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id);
1450	<
1451	<	// While we still have unreported hits...
1452	<	while(curr_left_obs_order != VMAXINT32 ||
1453	<	curr_right_obs_order != VMAXINT32)
1454	<	{
1455	<	// Chew up left singletons
1456	<	while (curr_left_obs_order < curr_right_obs_order &&
1457	<	curr_left_obs_order != VMAXINT32)
1458	<	{
1459	<	HitsForRead best_hits;
1460	<	best_hits.insert_id = curr_left_obs_order;
1461	<	FragmentAlignmentGrade grade;
1462	<
1463	<	// Process hits for left singleton, select best alignments
1464	<	read_best_alignments(curr_left_hit_group, grade, best_hits, gtf_junctions);
1465	<	update_junctions(best_hits, junctions);
1436	>
1437	>	// events include splice junction, indels, and fusion points.
1438	>	struct ConsensusEventsWorker
1439	>	{
1440	>	void operator()()
1441	>	{
1442	>	ReadTable it;
1443	>	vector<BAMHitFactory*> hit_factories;
1444	>	hit_factories.push_back(new BAMHitFactory(it, *rt));
1445	>	HitStream l_hs(left_map_fname, hit_factories.back(), false, true, true, true);
1446	>	if (left_map_offset > 0)
1447	>	l_hs.seek(left_map_offset);
1448	>
1449	>	hit_factories.push_back(new BAMHitFactory(it, *rt));
1450	>	HitStream r_hs(right_map_fname, hit_factories.back(), false, true, true, true);
1451	>	if (right_map_offset > 0)
1452	>	r_hs.seek(right_map_offset);
1453	>
1454	>	HitsForRead curr_left_hit_group;
1455	>	HitsForRead curr_right_hit_group;
1456	>
1457	>	l_hs.next_read_hits(curr_left_hit_group);
1458	>	r_hs.next_read_hits(curr_right_hit_group);
1459	>
1460	>	uint32_t curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id);
1461	>	uint32_t curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id);
1462	>
1463	>	// While we still have unreported hits...
1464	>	while((curr_left_obs_order != VMAXINT32 || curr_right_obs_order != VMAXINT32) &&
1465	>	(curr_left_obs_order < end_id || curr_right_obs_order < end_id))
1466	>	{
1467	>	// Chew up left singletons
1468	>	while (curr_left_obs_order < curr_right_obs_order &&
1469	>	curr_left_obs_order < end_id &&
1470	>	curr_left_obs_order != VMAXINT32)
1471	>	{
1472	>	HitsForRead best_hits;
1473	>	best_hits.insert_id = curr_left_obs_order;
1474
1475	<	// Get next hit group
1476	<	left_hs.next_read_hits(curr_left_hit_group);
1477	<	curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id);
1478	<	}
1475	>	// Process hits for left singleton, select best alignments
1476	>	read_best_alignments(curr_left_hit_group, best_hits, *gtf_junctions);
1477	>	update_coverage(best_hits, *coverage);
1478	>	update_junctions(best_hits, *junctions);
1479	>	update_insertions_and_deletions(best_hits, *insertions, *deletions);
1480	>	update_fusions(best_hits, *rt, *fusions);
1481	>
1482	>	// Get next hit group
1483	>	l_hs.next_read_hits(curr_left_hit_group);
1484	>	curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id);
1485	>	}
1486
1487	<	// Chew up right singletons
1488	<	while (curr_left_obs_order > curr_right_obs_order &&
1489	<	curr_right_obs_order != VMAXINT32)
1490	<	{
1491	<	HitsForRead best_hits;
1492	<	best_hits.insert_id = curr_right_obs_order;
1493	<	FragmentAlignmentGrade grade;
1494	<
1495	<	// Process hit for right singleton, select best alignments
1496	<	read_best_alignments(curr_right_hit_group,grade, best_hits, gtf_junctions);
1497	<	update_junctions(best_hits, junctions);
1487	>	// Chew up right singletons
1488	>	while (curr_left_obs_order > curr_right_obs_order &&
1489	>	curr_right_obs_order < end_id &&
1490	>	curr_right_obs_order != VMAXINT32)
1491	>	{
1492	>	HitsForRead best_hits;
1493	>	best_hits.insert_id = curr_right_obs_order;
1494	>	if (curr_right_obs_order >= begin_id)
1495	>	{
1496	>	// Process hit for right singleton, select best alignments
1497	>	read_best_alignments(curr_right_hit_group, best_hits, *gtf_junctions);
1498	>	update_coverage(best_hits, *coverage);
1499	>	update_junctions(best_hits, *junctions);
1500	>	update_insertions_and_deletions(best_hits, *insertions, *deletions);
1501	>	update_fusions(best_hits, *rt, *fusions);
1502	>	}
1503
1504	<	// Get next hit group
1505	<	right_hs.next_read_hits(curr_right_hit_group);
1506	<	curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id);
1507	<	}
1504	>	// Get next hit group
1505	>	r_hs.next_read_hits(curr_right_hit_group);
1506	>	curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id);
1507	>	}
1508
1509	<	// Since we have both left hits and right hits for this insert,
1510	<	// Find the best pairing and print both
1511	<	while (curr_left_obs_order == curr_right_obs_order &&
1512	<	curr_left_obs_order != VMAXINT32 && curr_right_obs_order != VMAXINT32)
1513	<	{
1514	<	if (curr_left_hit_group.hits.empty())
1515	<	{
1516	<	HitsForRead right_best_hits;
1517	<	right_best_hits.insert_id = curr_right_obs_order;
1518	<
1519	<	FragmentAlignmentGrade grade;
1520	<	read_best_alignments(curr_right_hit_group, grade, right_best_hits, gtf_junctions);
1521	<	update_junctions(right_best_hits, junctions);
1522	<	}
1523	<	else if (curr_right_hit_group.hits.empty())
1524	<	{
1525	<	HitsForRead left_best_hits;
1526	<	left_best_hits.insert_id = curr_left_obs_order;
1527	<
1528	<	FragmentAlignmentGrade grade;
1529	<	// Process hits for left singleton, select best alignments
1530	<	read_best_alignments(curr_left_hit_group, grade, left_best_hits, gtf_junctions);
1531	<	update_junctions(left_best_hits, junctions);
1532	<	}
1533	<	else
1534	<	{
1535	<	HitsForRead left_best_hits;
1536	<	HitsForRead right_best_hits;
1537	<	left_best_hits.insert_id = curr_left_obs_order;
1538	<	right_best_hits.insert_id = curr_right_obs_order;
1539	<
1540	<	// daehwan - apply gtf_junctions here, too!
1541	<
1542	<	InsertAlignmentGrade grade;
1543	<	pair_best_alignments(curr_left_hit_group,
1544	<	curr_right_hit_group,
1545	<	grade,
1546	<	left_best_hits,
1547	<	right_best_hits);
1548	<
1549	<	update_junctions(left_best_hits, junctions);
1550	<	update_junctions(right_best_hits, junctions);
1551	<	}
1552	<
705	<	left_hs.next_read_hits(curr_left_hit_group);
706	<	curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id);
1509	>	// Since we have both left hits and right hits for this insert,
1510	>	// Find the best pairing and print both
1511	>	while (curr_left_obs_order == curr_right_obs_order &&
1512	>	curr_left_obs_order < end_id &&
1513	>	curr_left_obs_order != VMAXINT32)
1514	>	{
1515	>	vector<pair<BowtieHit, BowtieHit> > best_hits;
1516	>
1517	>	InsertAlignmentGrade grade;
1518	>	pair_best_alignments(curr_left_hit_group,
1519	>	curr_right_hit_group,
1520	>	grade,
1521	>	best_hits,
1522	>	*gtf_junctions);
1523	>
1524	>	HitsForRead best_left_hit_group; best_left_hit_group.insert_id = curr_left_obs_order;
1525	>	HitsForRead best_right_hit_group; best_right_hit_group.insert_id = curr_left_obs_order;
1526	>
1527	>	if (best_hits.size() > 0)
1528	>	{
1529	>	for (size_t i = 0; i < best_hits.size(); ++i)
1530	>	{
1531	>	best_left_hit_group.hits.push_back(best_hits[i].first);
1532	>	best_right_hit_group.hits.push_back(best_hits[i].second);
1533	>	}
1534	>	}
1535	>	else
1536	>	{
1537	>	best_left_hit_group.hits = curr_left_hit_group.hits;
1538	>	best_right_hit_group.hits = curr_right_hit_group.hits;
1539	>	}
1540	>
1541	>	update_coverage(best_left_hit_group, *coverage);
1542	>	update_junctions(best_left_hit_group, *junctions);
1543	>	update_insertions_and_deletions(best_left_hit_group, *insertions, *deletions);
1544	>	update_fusions(best_left_hit_group, *rt, *fusions);
1545	>
1546	>	update_coverage(best_right_hit_group, *coverage);
1547	>	update_junctions(best_right_hit_group, *junctions);
1548	>	update_insertions_and_deletions(best_right_hit_group, *insertions, *deletions);
1549	>	update_fusions(best_right_hit_group, *rt, *fusions);
1550	>
1551	>	l_hs.next_read_hits(curr_left_hit_group);
1552	>	curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id);
1553
1554	<	right_hs.next_read_hits(curr_right_hit_group);
1555	<	curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id);
1556	<	}
1557	<	}
1558	<	left_hs.reset();
1559	<	right_hs.reset();
1560	<	}
1554	>	r_hs.next_read_hits(curr_right_hit_group);
1555	>	curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id);
1556	>	}
1557	>	}
1558	>
1559	>	for (size_t i = 0; i < hit_factories.size(); ++i)
1560	>	delete hit_factories[i];
1561	>
1562	>	hit_factories.clear();
1563	>	}
1564	>
1565	>	string left_map_fname;
1566	>	string right_map_fname;
1567	>
1568	>	RefSequenceTable* rt;
1569	>
1570	>	JunctionSet* gtf_junctions;
1571	>
1572	>	uint64_t begin_id;
1573	>	uint64_t end_id;
1574	>	int64_t left_map_offset;
1575	>	int64_t right_map_offset;
1576	>
1577	>	JunctionSet* junctions;
1578	>	InsertionSet* insertions;
1579	>	DeletionSet* deletions;
1580	>	FusionSet* fusions;
1581	>
1582	>	Coverage* coverage;
1583	>	};
1584	>
1585	>	struct ReportWorker
1586	>	{
1587	>	void write_singleton_alignments(uint64_t curr_obs_order,
1588	>	HitsForRead& curr_hit_group,
1589	>	ReadStream& reads_file,
1590	>	GBamWriter& bam_writer,
1591	>	GBamWriter& um_out,
1592	>	FragmentType fragment_type)
1593	>	{
1594	>	HitsForRead best_hits;
1595	>	best_hits.insert_id = curr_obs_order;
1596	>
1597	>	realign_reads(curr_hit_group, *rt, *junctions, *rev_junctions,
1598	>	*insertions, *deletions, *rev_deletions, *fusions);
1599	>
1600	>	exclude_hits_on_filtered_junctions(*junctions, curr_hit_group);
1601	>
1602	>	// Process hits for singleton, select best alignments
1603	>	const bool final_report = true;
1604	>	read_best_alignments(curr_hit_group, best_hits, *gtf_junctions,
1605	>	*junctions, *insertions, *deletions, *fusions, *coverage,
1606	>	final_report);
1607	>
1608	>	if (best_hits.hits.size() > 0)
1609	>	{
1610	>	Read read;
1611	>	bool got_read = reads_file.getRead(curr_obs_order, read,
1612	>	reads_format, false, begin_id, end_id,
1613	>	&um_out, false);
1614	>	assert (got_read);
1615	>	if (got_read)
1616	>	{
1617	>	update_junctions(best_hits, *final_junctions);
1618	>	update_insertions_and_deletions(best_hits, *final_insertions, *final_deletions);
1619	>	update_fusions(best_hits, *rt, *final_fusions, *fusions);
1620	>
1621	>	print_sam_for_single(*rt,
1622	>	best_hits,
1623	>	fragment_type,
1624	>	read.alt_name,
1625	>	bam_writer,
1626	>	rng);
1627	>	}
1628	>	}
1629	>	}
1630	>
1631	>	void operator()()
1632	>	{
1633	>	rng.seed(1);
1634	>
1635	>	ReadTable it;
1636	>	GBamWriter bam_writer(bam_output_fname.c_str(), sam_header.c_str());
1637	>
1638	>	ReadStream left_reads_file(left_reads_fname);
1639	>	if (left_reads_file.file() == NULL)
1640	>	err_die("Error: cannot open %s for reading\n", left_reads_fname.c_str());
1641	>	if (left_reads_file.isBam()) {
1642	>	left_reads_file.use_alt_name();
1643	>	left_reads_file.ignoreQC();
1644	>	}
1645	>	if (left_reads_offset > 0)
1646	>	left_reads_file.seek(left_reads_offset);
1647	>
1648	>	//if (!zpacker.empty()) left_um_fname+=".z";
1649	>	GBamWriter* left_um_out=new GBamWriter(left_um_fname.c_str(), sam_header.c_str());
1650	>	GBamWriter* right_um_out=NULL;
1651	>
1652	>	ReadStream right_reads_file(right_reads_fname);
1653	>	if (right_reads_offset > 0)
1654	>	right_reads_file.seek(right_reads_offset);
1655	>
1656	>	//FZPipe right_um_out;
1657	>	if (!right_reads_fname.empty())
1658	>	{
1659	>	if (right_reads_file.isBam()) {
1660	>	right_reads_file.use_alt_name();
1661	>	right_reads_file.ignoreQC();
1662	>	right_um_out=new GBamWriter(right_um_fname.c_str(), sam_header.c_str());
1663	>	}
1664	>	//if (!zpacker.empty()) right_um_fname+=".z";
1665	>	//if (right_um_out.openWrite(right_um_fname.c_str(), zpacker)==NULL)
1666	>	//  err_die("Error: cannot open file %s for writing!\n",right_um_fname.c_str());
1667	>	}
1668	>
1669	>	vector<BAMHitFactory*> hit_factories;
1670	>	hit_factories.push_back(new BAMHitFactory(it, *rt));
1671	>	HitStream left_hs(left_map_fname, hit_factories.back(), false, true, true, true);
1672	>	if (left_map_offset > 0)
1673	>	left_hs.seek(left_map_offset);
1674	>
1675	>	hit_factories.push_back(new BAMHitFactory(it, *rt));
1676	>	HitStream right_hs(right_map_fname, hit_factories.back(), false, true, true, true);
1677	>	if (right_map_offset > 0)
1678	>	right_hs.seek(right_map_offset);
1679	>
1680	>	HitsForRead curr_left_hit_group;
1681	>	HitsForRead curr_right_hit_group;
1682	>
1683	>	left_hs.next_read_hits(curr_left_hit_group);
1684	>	right_hs.next_read_hits(curr_right_hit_group);
1685	>
1686	>	uint64_t curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id);
1687	>	uint64_t curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id);
1688	>
1689	>	const bool final_report = true;
1690	>
1691	>	// While we still have unreported hits...
1692	>	Read l_read;
1693	>	Read r_read;
1694	>	while((curr_left_obs_order != VMAXINT32 || curr_right_obs_order != VMAXINT32) &&
1695	>	(curr_left_obs_order < end_id || curr_right_obs_order < end_id))
1696	>	{
1697	>	// Chew up left singletons (pairs with right reads unmapped)
1698	>	while (curr_left_obs_order < curr_right_obs_order &&
1699	>	curr_left_obs_order < end_id &&
1700	>	curr_left_obs_order != VMAXINT32)
1701	>	{
1702	>	write_singleton_alignments(curr_left_obs_order,
1703	>	curr_left_hit_group,
1704	>	left_reads_file,
1705	>	bam_writer,
1706	>	*left_um_out,
1707	>	right_map_fname.empty() ? FRAG_UNPAIRED : FRAG_LEFT);
1708	>
1709	>	// Get next hit group
1710	>	left_hs.next_read_hits(curr_left_hit_group);
1711	>	curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id);
1712	>	} //left singletons
1713	>
1714	>	// Chew up right singletons
1715	>	while (curr_left_obs_order > curr_right_obs_order &&
1716	>	curr_right_obs_order < end_id &&
1717	>	curr_right_obs_order != VMAXINT32)
1718	>	{
1719	>	write_singleton_alignments(curr_right_obs_order,
1720	>	curr_right_hit_group,
1721	>	right_reads_file,
1722	>	bam_writer,
1723	>	*right_um_out,
1724	>	FRAG_RIGHT);
1725	>
1726	>	// Get next hit group
1727	>	right_hs.next_read_hits(curr_right_hit_group);
1728	>	curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id);
1729	>	}
1730	>
1731	>	// Since we have both left hits and right hits for this insert,
1732	>	// Find the best pairing and print both
1733	>	while (curr_left_obs_order == curr_right_obs_order &&
1734	>	curr_left_obs_order < end_id &&
1735	>	curr_left_obs_order != VMAXINT32)
1736	>	{
1737	>	realign_reads(curr_left_hit_group, *rt, *junctions, *rev_junctions,
1738	>	*insertions, *deletions, *rev_deletions, *fusions);
1739	>	exclude_hits_on_filtered_junctions(*junctions, curr_left_hit_group);
1740	>
1741	>	realign_reads(curr_right_hit_group, *rt, *junctions, *rev_junctions,
1742	>	*insertions, *deletions, *rev_deletions, *fusions);
1743	>	exclude_hits_on_filtered_junctions(*junctions, curr_right_hit_group);
1744	>
1745	>	vector<pair<BowtieHit, BowtieHit> > best_hits;
1746
1747	+	bool paired_alignments = curr_left_hit_group.hits.size() > 0 && curr_right_hit_group.hits.size() > 0;
1748	+	InsertAlignmentGrade grade;
1749
1750	<	void driver(GBamWriter& bam_writer,
1751	<	string& left_map_fname,
1752	<	FILE* left_reads,
1753	<	string& right_map_fname,
1754	<	FILE* right_reads,
1750	>	if (paired_alignments)
1751	>	{
1752	>	pair_best_alignments(curr_left_hit_group,
1753	>	curr_right_hit_group,
1754	>	grade, best_hits, *gtf_junctions,
1755	>	*junctions, *insertions, *deletions, *fusions, *coverage,
1756	>	final_report);
1757	>
1758	>	if (best_hits.size() <= 0 ||
1759	>	(grade.fusion && !fusion_search && !report_discordant_pair_alignments))
1760	>	paired_alignments = false;
1761	>	}
1762	>
1763	>	if (paired_alignments)
1764	>	{
1765	>	HitsForRead best_left_hit_group; best_left_hit_group.insert_id = curr_left_obs_order;
1766	>	HitsForRead best_right_hit_group; best_right_hit_group.insert_id = curr_left_obs_order;
1767	>
1768	>	for (size_t i = 0; i < best_hits.size(); ++i)
1769	>	{
1770	>	best_left_hit_group.hits.push_back(best_hits[i].first);
1771	>	best_right_hit_group.hits.push_back(best_hits[i].second);
1772	>	}
1773	>
1774	>	if (best_hits.size() > 0)
1775	>	{
1776	>	bool got_left_read = left_reads_file.getRead(best_hits[0].first.insert_id(), l_read,
1777	>	reads_format, false, begin_id, end_id,
1778	>	left_um_out, false);
1779	>
1780	>	bool got_right_read = right_reads_file.getRead(best_hits[0].first.insert_id(), r_read,
1781	>	reads_format, false, begin_id, end_id,
1782	>	right_um_out, false);
1783	>
1784	>	assert (got_left_read && got_right_read);
1785	>
1786	>	if (got_left_read && got_right_read)
1787	>	{
1788	>	update_junctions(best_left_hit_group, *final_junctions);
1789	>	update_insertions_and_deletions(best_left_hit_group, *final_insertions, *final_deletions);
1790	>	update_fusions(best_left_hit_group, *rt, *final_fusions, *fusions);
1791	>
1792	>	update_junctions(best_right_hit_group, *final_junctions);
1793	>	update_insertions_and_deletions(best_right_hit_group, *final_insertions, *final_deletions);
1794	>	update_fusions(best_right_hit_group, *rt, *final_fusions, *fusions);
1795	>
1796	>	pair_support(best_hits, *final_fusions, *fusions);
1797	>
1798	>	print_sam_for_pair(*rt,
1799	>	best_hits,
1800	>	grade,
1801	>	bam_writer,
1802	>	l_read.alt_name,
1803	>	r_read.alt_name,
1804	>	rng,
1805	>	begin_id,
1806	>	end_id);
1807	>	}
1808	>	}
1809	>	}
1810	>	else
1811	>	{
1812	>	if (curr_left_hit_group.hits.size() > 0)
1813	>	{
1814	>	write_singleton_alignments(curr_left_obs_order,
1815	>	curr_left_hit_group,
1816	>	left_reads_file,
1817	>	bam_writer,
1818	>	*left_um_out,
1819	>	(right_map_fname.empty() ? FRAG_UNPAIRED : FRAG_LEFT));
1820	>	}
1821	>
1822	>	if (curr_right_hit_group.hits.size() > 0)
1823	>	{   //only right read mapped
1824	>	//write it in the mapped file with the #MAPPED# flag
1825	>	write_singleton_alignments(curr_right_obs_order,
1826	>	curr_right_hit_group,
1827	>	right_reads_file,
1828	>	bam_writer,
1829	>	*right_um_out,
1830	>	FRAG_RIGHT);
1831	>	}
1832	>	}
1833	>
1834	>	left_hs.next_read_hits(curr_left_hit_group);
1835	>	curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id);
1836	>
1837	>	right_hs.next_read_hits(curr_right_hit_group);
1838	>	curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id);
1839	>	}
1840	>	} //while we still have unreported hits..
1841	>	//print the remaining unmapped reads at the end of each reads' stream
1842	>
1843	>	left_reads_file.getRead(VMAXINT32, l_read,
1844	>	reads_format, false, begin_id, end_id,
1845	>	left_um_out, false);
1846	>	if (right_reads_file.file())
1847	>	right_reads_file.getRead(VMAXINT32, r_read,
1848	>	reads_format, false, begin_id, end_id,
1849	>	right_um_out, false);
1850	>
1851	>
1852	>	// pclose (pipe close), which waits for a process to end, seems to conflict with boost::thread::join somehow,
1853	>	// resulting in deadlock like behavior.
1854	>	delete left_um_out;
1855	>	delete right_um_out;
1856	>
1857	>
1858	>	for (size_t i = 0; i < hit_factories.size(); ++i)
1859	>	delete hit_factories[i];
1860	>
1861	>	hit_factories.clear();
1862	>	}
1863	>
1864	>	string bam_output_fname;
1865	>	string sam_header_fname;
1866	>
1867	>	string left_reads_fname;
1868	>	string left_map_fname;
1869	>	string right_reads_fname;
1870	>	string right_map_fname;
1871	>
1872	>	string left_um_fname;
1873	>	string right_um_fname;
1874	>
1875	>	JunctionSet* gtf_junctions;
1876	>
1877	>	JunctionSet* junctions;
1878	>	JunctionSet* rev_junctions;
1879	>	InsertionSet* insertions;
1880	>	DeletionSet* deletions;
1881	>	DeletionSet* rev_deletions;
1882	>	FusionSet* fusions;
1883	>	Coverage* coverage;
1884	>
1885	>	JunctionSet* final_junctions;
1886	>	InsertionSet* final_insertions;
1887	>	DeletionSet* final_deletions;
1888	>	FusionSet* final_fusions;
1889	>
1890	>	RefSequenceTable* rt;
1891	>
1892	>	uint64_t begin_id;
1893	>	uint64_t end_id;
1894	>	int64_t left_reads_offset;
1895	>	int64_t left_map_offset;
1896	>	int64_t right_reads_offset;
1897	>	int64_t right_map_offset;
1898	>
1899	>	boost::random::mt19937 rng;
1900	>	};
1901	>
1902	>	void driver(const string& bam_output_fname,
1903	>	istream& ref_stream,
1904	>	const string& left_map_fname,
1905	>	const string& left_reads_fname,
1906	>	const string& right_map_fname,
1907	>	const string& right_reads_fname,
1908		FILE* junctions_out,
1909		FILE* insertions_out,
1910		FILE* deletions_out,
1911	<	FILE* left_um_out,
726	<	FILE* right_um_out
727	<	)
1911	>	FILE* fusions_out)
1912		{
1913	<	ReadTable it;
1913	>	if (!parallel)
1914	>	num_threads = 1;
1915	>
1916		RefSequenceTable rt(sam_header, true);
1917	<	srandom(1);
1917	>	get_seqs(ref_stream, rt, true);
1918	>
1919	>	srandom(1);
1920	>
1921		JunctionSet gtf_junctions;
1922		if (!gtf_juncs.empty())
1923		{
1924	<	char splice_buf[2048];
1924	>	char splice_buf[4096];
1925		FILE* splice_coords = fopen(gtf_juncs.c_str(), "r");
1926		if (splice_coords)
1927		{
#	Line 759 | Line 1948
1948		uint32_t right_coord = atoi(scan_right_coord);
1949		bool antisense = *orientation == '-';
1950
1951	<	// add 1 to left_coord to meet BED format
1952	<	gtf_junctions.insert(make_pair<Junction, JunctionStats>(Junction(ref_id, left_coord + 1, right_coord, antisense), JunctionStats()));
1951	>	JunctionStats junction_stat;
1952	>	junction_stat.gtf_match = true;
1953	>	junction_stat.accepted = true;
1954	>
1955	>	gtf_junctions.insert(make_pair<Junction, JunctionStats>(Junction(ref_id, left_coord, right_coord, antisense), junction_stat));
1956		}
1957		}
1958		fprintf(stderr, "Loaded %d GFF junctions from %s.\n", (int)(gtf_junctions.size()), gtf_juncs.c_str());
1959		}
1960
1961	<	BAMHitFactory hit_factory(it,rt);
1962	<	JunctionSet junctions;
1963	<	{
772	<	HitStream l_hs(left_map_fname, &hit_factory, false, true, true, true);
773	<	HitStream r_hs(right_map_fname, &hit_factory, false, true, true, true);
774	<	get_junctions_from_best_hits(l_hs, r_hs, it, junctions, gtf_junctions);
775	<	//this resets the streams
776	<	}
777	<
778	<	HitStream left_hs(left_map_fname, &hit_factory, false, true, true, true);
779	<	HitStream right_hs(right_map_fname, &hit_factory, false, true, true, true);
780	<
781	<	size_t num_unfiltered_juncs = junctions.size();
782	<	fprintf(stderr, "Loaded %lu junctions\n", (long unsigned int) num_unfiltered_juncs);
783	<
784	<	HitsForRead curr_left_hit_group;
785	<	HitsForRead curr_right_hit_group;
786	<
787	<	left_hs.next_read_hits(curr_left_hit_group);
788	<	right_hs.next_read_hits(curr_right_hit_group);
789	<
790	<	uint32_t curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id);
791	<	uint32_t curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id);
792	<
793	<	// Read hits, extract junctions, and toss the ones that arent strongly enough supported.
794	<	filter_junctions(junctions, gtf_junctions);
795	<	//size_t num_juncs_after_filter = junctions.size();
796	<	//fprintf(stderr, "Filtered %lu junctions\n",
797	<	//     num_unfiltered_juncs - num_juncs_after_filter);
798	<
799	<	/*
800	<	size_t small_overhangs = 0;
801	<	for (JunctionSet::iterator i = junctions.begin(); i != junctions.end(); ++i)
1961	>	vector<uint64_t> read_ids;
1962	>	vector<vector<int64_t> > offsets;
1963	>	if (num_threads > 1)
1964		{
1965	<	if (i->second.accepted &&
1966	<	(i->second.left_extent < min_anchor_len || i->second.right_extent < min_anchor_len))
1967	<	{
1968	<	small_overhangs++;
1969	<	}
1965	>	vector<string> fnames;
1966	>	if (right_map_fname != "")
1967	>	{
1968	>	fnames.push_back(right_reads_fname);
1969	>	fnames.push_back(right_map_fname);
1970	>	}
1971	>	fnames.push_back(left_reads_fname);
1972	>	fnames.push_back(left_map_fname);
1973	>	bool enough_data = calculate_offsets(fnames, read_ids, offsets);
1974	>	if (!enough_data)
1975	>	num_threads = 1;
1976		}
809	–
810	–	if (small_overhangs >0)
811	–	fprintf(stderr, "Warning: %lu small overhang junctions!\n", (long unsigned int)small_overhangs);
812	–	*/
813	–	JunctionSet final_junctions; // the junctions formed from best hits
814	–	InsertionSet final_insertions;
815	–	DeletionSet final_deletions;
816	–
817	–	fprintf (stderr, "Reporting final accepted alignments...");
818	–	// While we still have unreported hits...
819	–	Read l_read;
820	–	Read r_read;
821	–	while(curr_left_obs_order != VMAXINT32 ||
822	–	curr_right_obs_order != VMAXINT32)
823	–	{
824	–	// Chew up left singletons (pairs with right reads unmapped)
825	–	while (curr_left_obs_order < curr_right_obs_order &&
826	–	curr_left_obs_order != VMAXINT32)
827	–	{
828	–	HitsForRead best_hits;
829	–	best_hits.insert_id = curr_left_obs_order;
830	–	FragmentAlignmentGrade grade;
831	–	bool got_read=get_read_from_stream(curr_left_obs_order,
832	–	left_reads,reads_format, false, l_read, left_um_out);
833	–	assert(got_read);
834	–	if (right_reads) {
835	–	fprintf(left_um_out, "@%s #MAPPED#\n%s\n+\n%s\n", l_read.alt_name.c_str(),
836	–	l_read.seq.c_str(), l_read.qual.c_str());
837	–	got_read=get_read_from_stream(curr_left_obs_order,
838	–	right_reads,reads_format, false,
839	–	r_read, right_um_out, true);
840	–	assert(got_read);
841	–	}
842	–	exclude_hits_on_filtered_junctions(junctions, curr_left_hit_group);
843	–
844	–	// Process hits for left singleton, select best alignments
845	–	read_best_alignments(curr_left_hit_group, grade, best_hits, gtf_junctions);
846	–	if (best_hits.hits.size()>0 && best_hits.hits.size() <= max_multihits)
847	–	{
848	–	update_junctions(best_hits, final_junctions);
849	–	update_insertions_and_deletions(best_hits, final_insertions, final_deletions);
850	–
851	–	print_sam_for_single(rt,
852	–	best_hits,
853	–	grade,
854	–	(right_map_fname.empty() ? FRAG_UNPAIRED : FRAG_LEFT),
855	–	l_read,
856	–	bam_writer,
857	–	left_um_out);
858	–	}
859	–
860	–	// Get next hit group
861	–	left_hs.next_read_hits(curr_left_hit_group);
862	–	curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id);
863	–	} //left singletons
864	–
865	–	// Chew up right singletons
866	–	while (curr_left_obs_order > curr_right_obs_order &&
867	–	curr_right_obs_order != VMAXINT32)
868	–	{
869	–	HitsForRead best_hits;
870	–	best_hits.insert_id = curr_right_obs_order;
871	–	FragmentAlignmentGrade grade;
872	–
873	–	bool got_read=get_read_from_stream(curr_right_obs_order,
874	–	right_reads,reads_format, false, r_read, right_um_out);
875	–	assert(got_read);
876	–	fprintf(right_um_out, "@%s #MAPPED#\n%s\n+\n%s\n", r_read.alt_name.c_str(),
877	–	r_read.seq.c_str(), r_read.qual.c_str());
878	–	got_read=get_read_from_stream(curr_right_obs_order,
879	–	left_reads,reads_format, false,
880	–	l_read, left_um_out, true);
881	–	assert(got_read);
1977
1978	<	exclude_hits_on_filtered_junctions(junctions, curr_right_hit_group);
1979	<
1980	<	// Process hit for right singleton, select best alignments
1981	<	read_best_alignments(curr_right_hit_group, grade, best_hits, gtf_junctions);
1982	<
1983	<	if (best_hits.hits.size()>0 && best_hits.hits.size() <= max_multihits)
1984	<	{
1985	<	update_junctions(best_hits, final_junctions);
1986	<	update_insertions_and_deletions(best_hits, final_insertions, final_deletions);
1987	<
893	<	print_sam_for_single(rt,
894	<	best_hits,
895	<	grade,
896	<	FRAG_RIGHT,
897	<	r_read,
898	<	bam_writer,
899	<	right_um_out);
900	<	}
901	<
902	<	// Get next hit group
903	<	right_hs.next_read_hits(curr_right_hit_group);
904	<	curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id);
905	<	}
906	<
907	<	// Since we have both left hits and right hits for this insert,
908	<	// Find the best pairing and print both
909	<	while (curr_left_obs_order == curr_right_obs_order &&
910	<	curr_left_obs_order != VMAXINT32 && curr_right_obs_order != VMAXINT32)
911	<	{
912	<	exclude_hits_on_filtered_junctions(junctions, curr_left_hit_group);
913	<	exclude_hits_on_filtered_junctions(junctions, curr_right_hit_group);
914	<
915	<	if (curr_left_hit_group.hits.empty())
916	<	{   //only right read mapped
917	<	//write it in the mapped file with the #MAPPED# flag
918	<
919	<	bool got_read=get_read_from_stream(curr_left_obs_order,
920	<	right_reads,reads_format, false, r_read, right_um_out);
921	<	assert(got_read);
922	<	fprintf(right_um_out, "@%s #MAPPED#\n%s\n+\n%s\n", r_read.alt_name.c_str(),
923	<	r_read.seq.c_str(), r_read.qual.c_str());
924	<	HitsForRead right_best_hits;
925	<	right_best_hits.insert_id = curr_right_obs_order;
926	<
927	<	FragmentAlignmentGrade grade;
928	<	read_best_alignments(curr_right_hit_group, grade, right_best_hits, gtf_junctions);
929	<
930	<	if (right_best_hits.hits.size()>0 && right_best_hits.hits.size() <= max_multihits)
931	<	{
932	<	update_junctions(right_best_hits, final_junctions);
933	<	update_insertions_and_deletions(right_best_hits, final_insertions, final_deletions);
934	<
935	<	print_sam_for_single(rt,
936	<	right_best_hits,
937	<	grade,
938	<	FRAG_RIGHT,
939	<	r_read,
940	<	bam_writer,
941	<	right_um_out);
942	<	}
943	<	}
944	<	else if (curr_right_hit_group.hits.empty())
945	<	{
946	<	HitsForRead left_best_hits;
947	<	left_best_hits.insert_id = curr_left_obs_order;
948	<	//only left read mapped
949	<	bool got_read=get_read_from_stream(curr_left_obs_order,
950	<	left_reads,reads_format, false, l_read, left_um_out);
951	<	assert(got_read);
952	<	fprintf(left_um_out, "@%s #MAPPED#\n%s\n+\n%s\n", l_read.alt_name.c_str(),
953	<	l_read.seq.c_str(), l_read.qual.c_str());
954	<	FragmentAlignmentGrade grade;
955	<	// Process hits for left singleton, select best alignments
956	<	read_best_alignments(curr_left_hit_group, grade, left_best_hits, gtf_junctions);
957	<
958	<	if (left_best_hits.hits.size()>0 && left_best_hits.hits.size() <= max_multihits)
959	<	{
960	<	update_junctions(left_best_hits, final_junctions);
961	<	update_insertions_and_deletions(left_best_hits, final_insertions, final_deletions);
962	<
963	<	print_sam_for_single(rt,
964	<	left_best_hits,
965	<	grade,
966	<	FRAG_LEFT,
967	<	l_read,
968	<	bam_writer,
969	<	left_um_out);
970	<	}
971	<	}
972	<	else
973	<	{   //hits for both left and right reads
974	<	HitsForRead left_best_hits;
975	<	HitsForRead right_best_hits;
976	<	left_best_hits.insert_id = curr_left_obs_order;
977	<	right_best_hits.insert_id = curr_right_obs_order;
978	<
979	<	InsertAlignmentGrade grade;
980	<	pair_best_alignments(curr_left_hit_group,
981	<	curr_right_hit_group,
982	<	grade,
983	<	left_best_hits,
984	<	right_best_hits);
985	<
986	<	if (left_best_hits.hits.size()>0 && left_best_hits.hits.size() <= max_multihits &&
987	<	right_best_hits.hits.size()>0 && right_best_hits.hits.size() <= max_multihits)
988	<	{
989	<	update_junctions(left_best_hits, final_junctions);
990	<	update_junctions(right_best_hits, final_junctions);
991	<	update_insertions_and_deletions(left_best_hits, final_insertions, final_deletions);
992	<	update_insertions_and_deletions(right_best_hits, final_insertions, final_deletions);
993	<
994	<	print_sam_for_pair(rt,
995	<	left_best_hits,
996	<	right_best_hits,
997	<	grade,
998	<	left_reads,
999	<	right_reads,
1000	<	bam_writer,
1001	<	left_um_out,
1002	<	right_um_out);
1003	<	}
1004	<	}
1005	<
1006	<	left_hs.next_read_hits(curr_left_hit_group);
1007	<	curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id);
1008	<
1009	<	right_hs.next_read_hits(curr_right_hit_group);
1010	<	curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id);
1011	<	}
1012	<
1013	<	} //while we still have unreported hits..
1014	<	//print the remaining unmapped reads at the end of each reads' stream
1015	<	get_read_from_stream(VMAXINT32,
1016	<	left_reads,
1017	<	reads_format,
1018	<	false,
1019	<	l_read,
1020	<	left_um_out);
1021	<	if (right_reads)
1022	<	get_read_from_stream(VMAXINT32,
1023	<	right_reads,
1024	<	reads_format,
1025	<	false,
1026	<	r_read,
1027	<	right_um_out);
1028	<	fprintf (stderr, "done.\n");
1029	<
1030	<	//small_overhangs = 0;
1031	<	for (JunctionSet::iterator i = final_junctions.begin(); i != final_junctions.end();)
1032	<	{
1033	<	if (i->second.supporting_hits == 0 || i->second.left_extent < 8 ||  i->second.right_extent < 8)
1034	<	{
1035	<	final_junctions.erase(i++);
1036	<	}
1037	<	else
1038	<	{
1039	<	++i;
1040	<	}
1041	<	}
1042	<
1043	<	//      if (small_overhangs > 0)
1044	<	//              fprintf(stderr, "Warning: %lu small overhang junctions!\n", small_overhangs);
1045	<
1046	<	fprintf (stderr, "Printing junction BED track...");
1047	<	print_junctions(junctions_out, final_junctions, rt);
1048	<	fprintf (stderr, "done\n");
1049	<
1050	<	fprintf (stderr, "Printing insertions...");
1051	<	print_insertions(insertions_out, final_insertions,rt);
1052	<	fclose(insertions_out);
1053	<	fprintf (stderr, "done\n");
1054	<
1055	<	fprintf (stderr, "Printing deletions...");
1056	<	print_deletions(deletions_out, final_deletions, rt);
1057	<	fclose(deletions_out);
1058	<	fprintf (stderr, "done\n");
1059	<
1060	<	fprintf(stderr, "Found %lu junctions from happy spliced reads\n", (long unsigned int)final_junctions.size());
1061	<	}
1062	<
1063	<	void print_usage()
1064	<	{
1065	<	fprintf(stderr, "Usage:   tophat_reports <junctions.bed> <insertions.vcf> <deletions.vcf> <accepted_hits.sam> <left_map1,...,left_mapN> <left_reads.fq>  [right_map1,...,right_mapN] [right_reads.fq]\n");
1066	<
1067	<	//      fprintf(stderr, "Usage:   tophat_reports <coverage.wig> <junctions.bed> <accepted_hits.sam> <map1.bwtout> [splice_map1.sbwtout]\n");
1068	<	}
1978	>	JunctionSet vjunctions[num_threads];
1979	>	InsertionSet vinsertions[num_threads];
1980	>	DeletionSet vdeletions[num_threads];
1981	>	FusionSet vfusions[num_threads];
1982	>	Coverage vcoverages[num_threads];
1983	>
1984	>	vector<boost::thread*> threads;
1985	>	for (int i = 0; i < num_threads; ++i)
1986	>	{
1987	>	ConsensusEventsWorker worker;
1988
1989	<	int main(int argc, char** argv)
1990	<	{
1991	<	fprintf(stderr, "tophat_reports v%s (%s)\n", PACKAGE_VERSION, SVN_REVISION);
1992	<	fprintf(stderr, "---------------------------------------\n");
1993	<
1994	<	reads_format = FASTQ;
1995	<
1996	<	int parse_ret = parse_options(argc, argv, print_usage);
1997	<	if (parse_ret)
1998	<	return parse_ret;
1999	<
2000	<	if(optind >= argc)
1989	>	worker.left_map_fname = left_map_fname;
1990	>	worker.right_map_fname = right_map_fname;
1991	>	worker.rt = &rt;
1992	>	worker.gtf_junctions = &gtf_junctions;
1993	>
1994	>	worker.junctions = &vjunctions[i];
1995	>	worker.insertions = &vinsertions[i];
1996	>	worker.deletions = &vdeletions[i];
1997	>	worker.fusions = &vfusions[i];
1998	>	worker.coverage = &vcoverages[i];
1999	>
2000	>	worker.right_map_offset = 0;
2001	>
2002	>	if (i == 0)
2003		{
2004	<	print_usage();
2005	<	return 1;
2004	>	worker.begin_id = 0;
2005	>	worker.left_map_offset = 0;
2006		}
2007	<
1087	<	string junctions_file_name = argv[optind++];
1088	<
1089	<	if(optind >= argc)
2007	>	else
2008		{
2009	<	print_usage();
2010	<	return 1;
2011	<	}
2012	<
2013	<	string insertions_file_name = argv[optind++];
2014	<	if(optind >= argc)
1097	<	{
1098	<	print_usage();
1099	<	return 1;
2009	>	size_t offsets_size = offsets[i-1].size();
2010	>
2011	>	worker.begin_id = read_ids[i-1];
2012	>	worker.left_map_offset = offsets[i-1].back();
2013	>	if (offsets_size == 4)
2014	>	worker.right_map_offset = offsets[i-1][1];
2015		}
2016	+
2017	+	worker.end_id = (i+1 < num_threads) ? read_ids[i] : std::numeric_limits<uint64_t>::max();
2018	+
2019	+	if (num_threads > 1)
2020	+	threads.push_back(new boost::thread(worker));
2021	+	else
2022	+	worker();
2023	+	}
2024	+
2025	+	for (size_t i = 0; i < threads.size(); ++i)
2026	+	{
2027	+	threads[i]->join();
2028	+	delete threads[i];
2029	+	threads[i] = NULL;
2030	+	}
2031	+	threads.clear();
2032	+
2033	+	JunctionSet& junctions = vjunctions[0];
2034	+	InsertionSet& insertions = vinsertions[0];
2035	+	DeletionSet& deletions = vdeletions[0];
2036	+	FusionSet& fusions = vfusions[0];
2037	+	Coverage& coverage = vcoverages[0];
2038	+	for (int i = 1; i < num_threads; ++i)
2039	+	{
2040	+	merge_with(junctions, vjunctions[i]);
2041	+	vjunctions[i].clear();
2042	+
2043	+	merge_with(insertions, vinsertions[i]);
2044	+	vinsertions[i].clear();
2045	+
2046	+	merge_with(deletions, vdeletions[i]);
2047	+	vdeletions[i].clear();
2048	+
2049	+	merge_with(fusions, vfusions[i]);
2050	+	vfusions[i].clear();
2051	+
2052	+	coverage.merge_with(vcoverages[i]);
2053	+	vcoverages[i].clear();
2054	+	}
2055	+
2056	+	merge_with(junctions, gtf_junctions);
2057	+
2058	+	coverage.calculate_coverage();
2059	+
2060	+	JunctionSet rev_junctions;
2061	+	JunctionSet::const_iterator junction_iter = junctions.begin();
2062	+	for (; junction_iter != junctions.end(); ++junction_iter)
2063	+	{
2064	+	const Junction& junction = junction_iter->first;
2065	+	Junction rev_junction = junction;
2066	+	rev_junction.left = junction.right;
2067	+	rev_junction.right = junction.left;
2068	+	rev_junctions[rev_junction] = junction_iter->second;
2069	+	}
2070	+
2071	+	DeletionSet rev_deletions;
2072	+	#if 0
2073	+	DeletionSet::const_iterator deletion_iter = deletions.begin();
2074	+	for (; deletion_iter != deletions.end(); ++deletion_iter)
2075	+	{
2076	+	const Deletion& deletion = deletion_iter->first;
2077	+	Deletion rev_deletion = deletion;
2078	+	rev_deletion.left = deletion.right;
2079	+	rev_deletion.right = deletion.left;
2080	+	rev_deletions[rev_deletion] = deletion_iter->second;
2081	+	}
2082	+	#endif
2083	+
2084	+	size_t num_unfiltered_juncs = junctions.size();
2085	+	fprintf(stderr, "Loaded %lu junctions\n", (long unsigned int) num_unfiltered_juncs);
2086	+
2087	+	// Read hits, extract junctions, and toss the ones that arent strongly enough supported.
2088	+	filter_junctions(junctions, gtf_junctions);
2089	+
2090	+	//size_t num_juncs_after_filter = junctions.size();
2091	+	//fprintf(stderr, "Filtered %lu junctions\n",
2092	+	//     num_unfiltered_juncs - num_juncs_after_filter);
2093	+
2094	+	/*
2095	+	size_t small_overhangs = 0;
2096	+	for (JunctionSet::iterator i = junctions.begin(); i != junctions.end(); ++i)
2097	+	{
2098	+	if (i->second.accepted &&
2099	+	(i->second.left_extent < min_anchor_len || i->second.right_extent < min_anchor_len))
2100	+	{
2101	+	small_overhangs++;
2102	+	}
2103	+	}
2104
2105	<	string deletions_file_name = argv[optind++];
2106	<	if(optind >= argc)
2105	>	if (small_overhangs >0)
2106	>	fprintf(stderr, "Warning: %lu small overhang junctions!\n", (long unsigned int)small_overhangs);
2107	>	*/
2108	>
2109	>	JunctionSet vfinal_junctions[num_threads];
2110	>	InsertionSet vfinal_insertions[num_threads];
2111	>	DeletionSet vfinal_deletions[num_threads];
2112	>	FusionSet vfinal_fusions[num_threads];
2113	>
2114	>	for (int i = 0; i < num_threads; ++i)
2115	>	{
2116	>	ReportWorker worker;
2117	>
2118	>	worker.sam_header_fname = sam_header;
2119	>	char filename[1024] = {0};
2120	>	sprintf(filename, "%s%d.bam", bam_output_fname.c_str(), i);
2121	>	worker.bam_output_fname = filename;
2122	>	string tmpoutdir=getFdir(worker.bam_output_fname);
2123	>	worker.left_um_fname = tmpoutdir;
2124	>	sprintf(filename, "unmapped_left_%d.bam", i);
2125	>	worker.left_um_fname+=filename;
2126	>
2127	>	if (right_reads_fname != "")
2128		{
2129	<	print_usage();
2130	<	return 1;
2129	>	sprintf(filename, "unmapped_right_%d.bam", i);
2130	>	worker.right_um_fname = tmpoutdir;
2131	>	worker.right_um_fname += filename;
2132		}
2133	<
2134	<	string accepted_hits_file_name = argv[optind++];
2135	<
2136	<	if(optind >= argc)
2133	>
2134	>	worker.left_reads_fname = left_reads_fname;
2135	>	worker.left_map_fname = left_map_fname;
2136	>	worker.right_reads_fname = right_reads_fname;
2137	>	worker.right_map_fname = right_map_fname;
2138	>
2139	>	worker.gtf_junctions = &gtf_junctions;
2140	>	worker.junctions = &junctions;
2141	>	worker.rev_junctions = &rev_junctions;
2142	>	worker.insertions = &insertions;
2143	>	worker.deletions = &deletions;
2144	>	worker.rev_deletions = &rev_deletions;
2145	>	worker.fusions = &fusions;
2146	>	worker.coverage = &coverage;
2147	>
2148	>	worker.final_junctions = &vfinal_junctions[i];
2149	>	worker.final_insertions = &vfinal_insertions[i];
2150	>	worker.final_deletions = &vfinal_deletions[i];
2151	>	worker.final_fusions = &vfinal_fusions[i];
2152	>	worker.rt = &rt;
2153	>
2154	>	worker.right_reads_offset = 0;
2155	>	worker.right_map_offset = 0;
2156	>
2157	>	if (i == 0)
2158		{
2159	<	print_usage();
2160	<	return 1;
2159	>	worker.begin_id = 0;
2160	>	worker.left_reads_offset = 0;
2161	>	worker.left_map_offset = 0;
2162		}
2163	<
1117	<	string left_map_filename = argv[optind++];
1118	<
1119	<	if(optind >= argc)
2163	>	else
2164		{
2165	<	print_usage();
2166	<	return 1;
2165	>	size_t offsets_size = offsets[i-1].size();
2166	>
2167	>	worker.begin_id = read_ids[i-1];
2168	>	worker.left_reads_offset = offsets[i-1][offsets_size - 2];
2169	>	worker.left_map_offset = offsets[i-1].back();
2170	>	if (offsets_size == 4)
2171	>	{
2172	>	worker.right_reads_offset = offsets[i-1][0];
2173	>	worker.right_map_offset = offsets[i-1][1];
2174	>	}
2175		}
2176	<	//FZPipe left_map_file;
2177	<	//string unbamcmd=getBam2SamCmd(left_map_filename);
1126	<	//left_map_file.openRead(left_map_filename, unbamcmd);
1127	<	string left_reads_filename = argv[optind++];
1128	<	string unzcmd=getUnpackCmd(left_reads_filename, false);
1129	<
1130	<	string right_map_filename;
1131	<	string right_reads_filename;
1132	<	FZPipe right_reads_file;
1133	<	string right_um_filename;
1134	<	FZPipe right_um_file;
2176	>
2177	>	worker.end_id = (i+1 < num_threads) ? read_ids[i] : std::numeric_limits<uint64_t>::max();
2178
2179	<	if (optind < argc)
2180	<	{
2181	<	right_map_filename = argv[optind++];
2182	<	if(optind >= argc) {
2183	<	print_usage();
1141	<	return 1;
1142	<	}
1143	<	//right_map_file.openRead(right_map_filename, unbamcmd);
1144	<	//if (optind<argc) {
1145	<	right_reads_filename=argv[optind++];
1146	<	right_reads_file.openRead(right_reads_filename,unzcmd);
1147	<	right_um_filename=output_dir+"/unmapped_right.fq";
1148	<	if (!zpacker.empty()) right_um_filename+=".z";
1149	<	if (right_um_file.openWrite(right_um_filename.c_str(), zpacker)==NULL)
1150	<	err_die("Error: cannot open file %s for writing!\n",right_um_filename.c_str());
2179	>	if (num_threads > 1)
2180	>	threads.push_back(new boost::thread(worker));
2181	>	else
2182	>	worker();
2183	>	}
2184
2185	<	//  }
2186	<	}
2187	<	FILE* junctions_file = fopen(junctions_file_name.c_str(), "w");
2188	<	if (junctions_file == NULL)
2189	<	{
2190	<	fprintf(stderr, "Error: cannot open BED file %s for writing\n",
2191	<	junctions_file_name.c_str());
2192	<	exit(1);
2193	<	}
2194	<
2195	<	FILE* insertions_file = fopen(insertions_file_name.c_str(), "w");
2196	<	if (insertions_file == NULL)
2185	>	for (size_t i = 0; i < threads.size(); ++i)
2186	>	{
2187	>	threads[i]->join();
2188	>	delete threads[i];
2189	>	threads[i] = NULL;
2190	>	}
2191	>	threads.clear();
2192	>
2193	>	JunctionSet final_junctions = vfinal_junctions[0];
2194	>	InsertionSet final_insertions = vfinal_insertions[0];
2195	>	DeletionSet final_deletions = vfinal_deletions[0];
2196	>	FusionSet final_fusions = vfinal_fusions[0];
2197	>	for (int i = 1; i < num_threads; ++i)
2198	>	{
2199	>	merge_with(final_junctions, vfinal_junctions[i]);
2200	>	vfinal_junctions[i].clear();
2201	>
2202	>	merge_with(final_insertions, vfinal_insertions[i]);
2203	>	vfinal_insertions[i].clear();
2204	>
2205	>	merge_with(final_deletions, vfinal_deletions[i]);
2206	>	vfinal_deletions[i].clear();
2207	>
2208	>	merge_with(final_fusions, vfinal_fusions[i]);
2209	>	vfinal_fusions[i].clear();
2210	>	}
2211	>
2212	>	fprintf (stderr, "Reporting final accepted alignments...");
2213	>	fprintf (stderr, "done.\n");
2214	>
2215	>	//small_overhangs = 0;
2216	>	for (JunctionSet::iterator i = final_junctions.begin(); i != final_junctions.end();)
2217	>	{
2218	>	if (i->second.supporting_hits == 0 || i->second.left_extent < 8 || i->second.right_extent < 8)
2219		{
2220	<	fprintf(stderr, "Error: cannot open VCF file %s for writing\n",
1166	<	insertions_file_name.c_str());
1167	<	exit(1);
2220	>	final_junctions.erase(i++);
2221		}
2222	<
1170	<
1171	<	FILE* deletions_file = fopen(deletions_file_name.c_str(), "w");
1172	<	if (deletions_file == NULL)
2222	>	else
2223		{
2224	<	fprintf(stderr, "Error: cannot open VCF file %s for writing\n",
1175	<	deletions_file_name.c_str());
1176	<	exit(1);
2224	>	++i;
2225		}
2226	<	bool uncompressed_bam=(accepted_hits_file_name=="-");
2227	<	GBamWriter bam_writer(accepted_hits_file_name.c_str(), sam_header.c_str(), uncompressed_bam);
2226	>	}
2227	>
2228	>	//    if (small_overhangs > 0)
2229	>	//            fprintf(stderr, "Warning: %lu small overhang junctions!\n", small_overhangs);
2230	>
2231	>	fprintf (stderr, "Printing junction BED track...");
2232	>	print_junctions(junctions_out, final_junctions, rt);
2233	>	fprintf (stderr, "done\n");
2234	>
2235	>	fprintf (stderr, "Printing insertions...");
2236	>	print_insertions(insertions_out, final_insertions,rt);
2237	>	fclose(insertions_out);
2238	>	fprintf (stderr, "done\n");
2239	>
2240	>	fprintf (stderr, "Printing deletions...");
2241	>	print_deletions(deletions_out, final_deletions, rt);
2242	>	fclose(deletions_out);
2243	>	fprintf (stderr, "done\n");
2244
2245	<	FZPipe left_reads_file(left_reads_filename, unzcmd);
2246	<	if (left_reads_file.file==NULL)
2247	<	{
2248	<	fprintf(stderr, "Error: cannot open reads file %s for reading\n",
2249	<	left_reads_filename.c_str());
2250	<	exit(1);
2251	<	}
2252	<	string left_um_filename(output_dir+"/unmapped_left.fq");
2253	<	if (!zpacker.empty()) left_um_filename+=".z";
1190	<	FZPipe left_um_file;
1191	<	if (left_um_file.openWrite(left_um_filename.c_str(), zpacker)==NULL)
1192	<	err_die("Error: cannot open file %s for writing!\n",left_um_filename.c_str());
1193	<
1194	<	driver(bam_writer, left_map_filename,
1195	<	left_reads_file.file,
1196	<	right_map_filename,
1197	<	right_reads_file.file,
1198	<	junctions_file,
1199	<	insertions_file,
1200	<	deletions_file,
1201	<	left_um_file.file,
1202	<	right_um_file.file);
1203	<	left_um_file.close();
1204	<	right_um_file.close();
1205	<	return 0;
2245	>	if (fusion_search)
2246	>	{
2247	>	fprintf (stderr, "Printing fusions...");
2248	>	print_fusions(fusions_out, final_fusions, rt);
2249	>	fclose(fusions_out);
2250	>	fprintf (stderr, "done\n");
2251	>	}
2252	>
2253	>	fprintf(stderr, "Found %lu junctions from happy spliced reads\n", (long unsigned int)final_junctions.size());
2254		}
2255
2256	+	void print_usage()
2257	+	{
2258	+	fprintf(stderr, "Usage:   tophat_reports <junctions.bed> <insertions.vcf> <deletions.vcf> <accepted_hits.sam> <left_map1,...,left_mapN> <left_reads.fq>  [right_map1,...,right_mapN] [right_reads.fq]\n");
2259	+
2260	+	//      fprintf(stderr, "Usage:   tophat_reports <coverage.wig> <junctions.bed> <accepted_hits.sam> <map1.bwtout> [splice_map1.sbwtout]\n");
2261	+	}
2262
2263	+	int main(int argc, char** argv)
2264	+	{
2265	+	fprintf(stderr, "tophat_reports v%s (%s)\n", PACKAGE_VERSION, SVN_REVISION);
2266	+	fprintf(stderr, "---------------------------------------\n");
2267	+
2268	+	reads_format = FASTQ;
2269	+
2270	+	int parse_ret = parse_options(argc, argv, print_usage);
2271	+	if (parse_ret)
2272	+	return parse_ret;
2273	+
2274	+	if(optind >= argc)
2275	+	{
2276	+	print_usage();
2277	+	return 1;
2278	+	}
2279	+
2280	+	string ref_file_name = argv[optind++];
2281	+
2282	+	if(optind >= argc)
2283	+	{
2284	+	print_usage();
2285	+	return 1;
2286	+	}
2287	+
2288	+	string junctions_file_name = argv[optind++];
2289	+
2290	+	if(optind >= argc)
2291	+	{
2292	+	print_usage();
2293	+	return 1;
2294	+	}
2295	+
2296	+	string insertions_file_name = argv[optind++];
2297	+	if(optind >= argc)
2298	+	{
2299	+	print_usage();
2300	+	return 1;
2301	+	}
2302	+
2303	+	string deletions_file_name = argv[optind++];
2304	+	if(optind >= argc)
2305	+	{
2306	+	print_usage();
2307	+	return 1;
2308	+	}
2309	+
2310	+	string fusions_file_name = argv[optind++];
2311	+	if(optind >= argc)
2312	+	{
2313	+	print_usage();
2314	+	return 1;
2315	+	}
2316	+
2317	+	string accepted_hits_file_name = argv[optind++];
2318	+	if(optind >= argc)
2319	+	{
2320	+	print_usage();
2321	+	return 1;
2322	+	}
2323	+
2324	+	string left_map_filename = argv[optind++];
2325	+	if(optind >= argc)
2326	+	{
2327	+	print_usage();
2328	+	return 1;
2329	+	}
2330	+
2331	+	string left_reads_filename = argv[optind++];
2332	+	string unzcmd=getUnpackCmd(left_reads_filename, false);
2333	+	string right_map_filename;
2334	+	string right_reads_filename;
2335	+
2336	+	if (optind < argc)
2337	+	{
2338	+	right_map_filename = argv[optind++];
2339	+	if(optind >= argc) {
2340	+	print_usage();
2341	+	return 1;
2342	+	}
2343	+	right_reads_filename=argv[optind++];
2344	+	}
2345	+
2346	+	ifstream ref_stream(ref_file_name.c_str(), ifstream::in);
2347	+	if (!ref_stream.good())
2348	+	{
2349	+	fprintf(stderr, "Error: cannot open %s for reading\n",
2350	+	ref_file_name.c_str());
2351	+	exit(1);
2352	+	}
2353	+
2354	+	FILE* junctions_file = fopen(junctions_file_name.c_str(), "w");
2355	+	if (junctions_file == NULL)
2356	+	{
2357	+	fprintf(stderr, "Error: cannot open BED file %s for writing\n",
2358	+	junctions_file_name.c_str());
2359	+	exit(1);
2360	+	}
2361	+
2362	+	FILE* insertions_file = fopen(insertions_file_name.c_str(), "w");
2363	+	if (insertions_file == NULL)
2364	+	{
2365	+	fprintf(stderr, "Error: cannot open VCF file %s for writing\n",
2366	+	insertions_file_name.c_str());
2367	+	exit(1);
2368	+	}
2369	+
2370	+	FILE* deletions_file = fopen(deletions_file_name.c_str(), "w");
2371	+	if (deletions_file == NULL)
2372	+	{
2373	+	fprintf(stderr, "Error: cannot open VCF file %s for writing\n",
2374	+	deletions_file_name.c_str());
2375	+	exit(1);
2376	+	}
2377	+
2378	+	FILE* fusions_file = NULL;
2379	+	if (fusion_search)
2380	+	{
2381	+	fusions_file = fopen(fusions_file_name.c_str(), "w");
2382	+	if (fusions_file == NULL)
2383	+	{
2384	+	fprintf(stderr, "Error: cannot open VCF file %s for writing\n",
2385	+	fusions_file_name.c_str());
2386	+	exit(1);
2387	+	}
2388	+	}
2389	+
2390	+	driver(accepted_hits_file_name,
2391	+	ref_stream,
2392	+	left_map_filename,
2393	+	left_reads_filename,
2394	+	right_map_filename,
2395	+	right_reads_filename,
2396	+	junctions_file,
2397	+	insertions_file,
2398	+	deletions_file,
2399	+	fusions_file);
2400	+
2401	+	return 0;
2402	+	}

Diff Legend

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