[ViewVC] Diff of: gclib/tophat_cpp/tophat

#	Line 73 \| Line 73
73		}
74		}
75
76	+	struct cmp_read_alignment
77	+	{
78	+	cmp_read_alignment(const JunctionSet& gtf_juncs) :
79	+	_gtf_juncs(&gtf_juncs) {}
80	+
81	+	bool operator() (const BowtieHit& l, const BowtieHit& r) const
82	+	{
83	+	FragmentAlignmentGrade gl(l, *_gtf_juncs);
84	+	FragmentAlignmentGrade gr(r, *_gtf_juncs);
85	+	bool better = gr < gl;
86	+	bool worse = gl < gr;
87	+
88	+	if (better && !worse)
89	+	return true;
90	+	else
91	+	return false;
92	+	}
93	+
94	+	const JunctionSet* _gtf_juncs;
95	+	};
96	+
97		void read_best_alignments(const HitsForRead& hits_for_read,
98	<	FragmentAlignmentGrade& best_grade,
99	<	HitsForRead& best_hits,
100	<	const JunctionSet& gtf_junctions)
98	>	FragmentAlignmentGrade& best_grade,
99	>	HitsForRead& best_hits,
100	>	const JunctionSet& gtf_junctions)
101		{
102		const vector<BowtieHit>& hits = hits_for_read.hits;
103		for (size_t i = 0; i < hits.size(); ++i)
104		{
105	<	if (hits[i].edit_dist()>max_read_mismatches) continue;
106	<
107	<	FragmentAlignmentGrade g(hits[i], gtf_junctions);
108	<
109	<	// Is the new status better than the current best one?
110	<	if (best_grade < g)
111	<	{
112	<	best_hits.hits.clear();
113	<	best_grade = g;
114	<	best_hits.hits.push_back(hits[i]);
115	<	}
116	<	else if (!(g < best_grade)) // is it just as good?
117	<	{
118	<	best_grade.num_alignments++;
119	<	best_hits.hits.push_back(hits[i]);
120	<	}
105	>	if (hits[i].edit_dist() > max_read_mismatches)
106	>	continue;
107	>
108	>	if (report_secondary_alignments)
109	>	{
110	>	best_hits.hits.push_back(hits[i]);
111	>	}
112	>	else
113	>	{
114	>	FragmentAlignmentGrade g(hits[i], gtf_junctions);
115	>
116	>	// Is the new status better than the current best one?
117	>	if (best_grade < g)
118	>	{
119	>	best_hits.hits.clear();
120	>	best_grade = g;
121	>	best_hits.hits.push_back(hits[i]);
122	>	}
123	>	else if (!(g < best_grade)) // is it just as good?
124	>	{
125	>	best_grade.num_alignments++;
126	>	best_hits.hits.push_back(hits[i]);
127	>	}
128	>	}
129	>	}
130	>
131	>	if (report_secondary_alignments && best_hits.hits.size() > 0)
132	>	{
133	>	cmp_read_alignment cmp(gtf_junctions);
134	>	sort(best_hits.hits.begin(), best_hits.hits.end(), cmp);
135	>	best_grade = FragmentAlignmentGrade(best_hits.hits[0], gtf_junctions);
136	>	best_grade.num_alignments = best_hits.hits.size();
137		}
138		}
139
140	<	void pair_best_alignments(const HitsForRead& left_hits,
104	<	const HitsForRead& right_hits,
105	<	InsertAlignmentGrade& best_grade,
106	<	HitsForRead& left_best_hits,
107	<	HitsForRead& right_best_hits)
140	>	void set_insert_alignment_grade(const BowtieHit& lh, const BowtieHit& rh, InsertAlignmentGrade& grade)
141		{
142		// max mate inner distance (genomic)
143		int min_mate_inner_dist = inner_dist_mean - inner_dist_std_dev;
144		if (max_mate_inner_dist == -1)
145		max_mate_inner_dist = inner_dist_mean + inner_dist_std_dev;
146
147	+	bool fusion = false;
148	+	if (fusion_search)
149	+	{
150	+	bool left_fusion = lh.fusion_opcode() != FUSION_NOTHING;
151	+	bool right_fusion = rh.fusion_opcode() != FUSION_NOTHING;
152	+	if (left_fusion && right_fusion)
153	+	return;
154	+
155	+	fusion = left_fusion \|\| right_fusion;
156	+	if (!fusion && lh.ref_id() != rh.ref_id())
157	+	fusion = true;
158	+
159	+	if (!fusion && lh.ref_id() == rh.ref_id())
160	+	{
161	+	if (lh.antisense_align() == rh.antisense_align())
162	+	fusion = true;
163	+	else
164	+	{
165	+	int inter_dist = 0;
166	+	if (lh.antisense_align())
167	+	inter_dist = lh.left() - rh.right();
168	+	else
169	+	inter_dist = rh.left() - lh.right();
170	+
171	+	if (inter_dist < -(int)max_insertion_length \|\| inter_dist > (int)fusion_min_dist)
172	+	fusion = true;
173	+	}
174	+	}
175	+	}
176	+
177	+	grade = InsertAlignmentGrade(lh, rh, min_mate_inner_dist, max_mate_inner_dist, fusion);
178	+	}
179	+
180	+	struct cmp_pair_alignment
181	+	{
182	+	cmp_pair_alignment(const JunctionSet& gtf_juncs) :
183	+	_gtf_juncs(&gtf_juncs) {}
184	+
185	+	bool operator() (const pair<BowtieHit, BowtieHit>& l, const pair<BowtieHit, BowtieHit>& r) const
186	+	{
187	+	InsertAlignmentGrade gl; set_insert_alignment_grade(l.first, l.second, gl);
188	+	InsertAlignmentGrade gr; set_insert_alignment_grade(r.first, r.second, gr);
189	+
190	+	bool better = gr < gl;
191	+	bool worse = gl < gr;
192	+
193	+	if (better && !worse)
194	+	return true;
195	+	else
196	+	return false;
197	+	}
198	+
199	+	const JunctionSet* _gtf_juncs;
200	+	};
201	+
202	+	void pair_best_alignments(const HitsForRead& left_hits,
203	+	const HitsForRead& right_hits,
204	+	InsertAlignmentGrade& best_grade,
205	+	vector<pair<BowtieHit, BowtieHit> >& best_hits,
206	+	const JunctionSet& gtf_junctions)
207	+	{
208		const vector<BowtieHit>& left = left_hits.hits;
209		const vector<BowtieHit>& right = right_hits.hits;
210
#	Line 123 \| Line 217
217		const BowtieHit& rh = right[j];
218		if (rh.edit_dist() > max_read_mismatches) continue;
219
220	<	bool fusion = false;
127	<	if (fusion_search)
128	<	{
129	<	bool left_fusion = lh.fusion_opcode() != FUSION_NOTHING;
130	<	bool right_fusion = rh.fusion_opcode() != FUSION_NOTHING;
131	<	if (left_fusion && right_fusion)
132	<	continue;
133	<
134	<	fusion = left_fusion \|\| right_fusion;
135	<	if (!fusion && lh.ref_id() != rh.ref_id())
136	<	fusion = true;
137	<
138	<	if (!fusion && lh.ref_id() == rh.ref_id())
139	<	{
140	<	if (lh.antisense_align() == rh.antisense_align())
141	<	fusion = true;
142	<	else
143	<	{
144	<	int inter_dist = 0;
145	<	if (lh.antisense_align())
146	<	inter_dist = lh.left() - rh.right();
147	<	else
148	<	inter_dist = rh.left() - lh.right();
149	<
150	<	if (inter_dist < -(int)max_insertion_length \|\| inter_dist > (int)fusion_min_dist)
151	<	fusion = true;
152	<	}
153	<	}
154	<	}
155	<
156	<	InsertAlignmentGrade g(lh, rh, min_mate_inner_dist, max_mate_inner_dist, fusion);
220	>	InsertAlignmentGrade g; set_insert_alignment_grade(lh, rh, g);
221
222	<	// Is the new status better than the current best one?
159	<	if (best_grade < g)
222	>	if (report_secondary_alignments)
223		{
224	<	left_best_hits.hits.clear();
162	<	right_best_hits.hits.clear();
163	<
164	<	best_grade = g;
165	<	left_best_hits.hits.push_back(lh);
166	<	right_best_hits.hits.push_back(rh);
224	>	best_hits.push_back(make_pair(lh, rh));
225		}
226	<	else if (!(g < best_grade))
226	>	else
227		{
228	<	left_best_hits.hits.push_back(lh);
229	<	right_best_hits.hits.push_back(rh);
230	<	best_grade.num_alignments++;
228	>	// Is the new status better than the current best one?
229	>	if (best_grade < g)
230	>	{
231	>	best_hits.clear();
232	>	best_grade = g;
233	>	best_hits.push_back(make_pair(lh, rh));
234	>	}
235	>	else if (!(g < best_grade))
236	>	{
237	>	best_hits.push_back(make_pair(lh, rh));
238	>	best_grade.num_alignments++;
239	>	}
240		}
241		}
242		}
243	+
244	+	if (report_secondary_alignments && best_hits.size() > 0)
245	+	{
246	+	cmp_pair_alignment cmp(gtf_junctions);
247	+	sort(best_hits.begin(), best_hits.end(), cmp);
248	+	set_insert_alignment_grade(best_hits[0].first, best_hits[0].second, best_grade);
249	+	best_grade.num_alignments = best_hits.size();
250	+	}
251		}
252
253		enum FragmentType {FRAG_UNPAIRED, FRAG_LEFT, FRAG_RIGHT};
#	Line 180 \| Line 255
255		void add_auxData(vector<string>& auxdata, vector<string>& sam_toks,
256		const RefSequenceTable& rt,const BowtieHit& bh, FragmentType insert_side,
257		int num_hits, const BowtieHit* next_hit, int hitIndex) {
258	+	bool XS_found = false;
259	+	if (sam_toks.size()>11) {
260
261	<	if (sam_toks.size()>11) {
262	<	for (size_t i=11;i<sam_toks.size();++i) {
263	<	if (sam_toks[i].find("NH:i:")==string::npos)
264	<	auxdata.push_back(sam_toks[i]);
265	<	}
266	<	}
267	<	string aux("NH:i:");
268	<	str_appendInt(aux, num_hits);
261	>	for (size_t i=11;i<sam_toks.size();++i) {
262	>	if (sam_toks[i].find("NH:i:")==string::npos &&
263	>	sam_toks[i].find("XS:i:")==string::npos)
264	>	auxdata.push_back(sam_toks[i]);
265	>
266	>	if (sam_toks[i].find("XS:A:")!=string::npos)
267	>	XS_found = true;
268	>	}
269	>	}
270	>	string aux("NH:i:");
271	>	str_appendInt(aux, num_hits);
272	>	auxdata.push_back(aux);
273	>	if (next_hit) {
274	>	const char* nh_ref_name = "=";
275	>	nh_ref_name = rt.get_name(next_hit->ref_id());
276	>	assert (nh_ref_name != NULL);
277	>	bool same_contig=(next_hit->ref_id()==bh.ref_id());
278	>	aux="CC:Z:";
279	>	aux+= (same_contig ? "=" : nh_ref_name);
280	>	auxdata.push_back(aux);
281	>	aux="CP:i:";
282	>	int nh_gpos=next_hit->left() + 1;
283	>	str_appendInt(aux, nh_gpos);
284		auxdata.push_back(aux);
285	<	if (next_hit) {
194	<	const char* nh_ref_name = "=";
195	<	nh_ref_name = rt.get_name(next_hit->ref_id());
196	<	assert (nh_ref_name != NULL);
197	<	bool same_contig=(next_hit->ref_id()==bh.ref_id());
198	<	aux="CC:Z:";
199	<	aux+= (same_contig ? "=" : nh_ref_name);
200	<	auxdata.push_back(aux);
201	<	aux="CP:i:";
202	<	int nh_gpos=next_hit->left() + 1;
203	<	str_appendInt(aux, nh_gpos);
204	<	auxdata.push_back(aux);
205	<	} //has next_hit
285	>	} //has next_hit
286		// FIXME: this code is still a bit brittle, because it contains no
287		// consistency check that the mates are on opposite strands - a current protocol
288		// requirement, and that the strand indicated by the alignment is consistent
289		// with the orientation of the splices (though that should be handled upstream).
290	+	if (!XS_found) {
291		const string xs_f("XS:A:+");
292		const string xs_r("XS:A:-");
293	<	if (bh.contiguous()) {
294	<	if (library_type == FR_FIRSTSTRAND) {
295	<	if (insert_side == FRAG_LEFT \|\| insert_side == FRAG_UNPAIRED) {
296	<	if (bh.antisense_align())
297	<	auxdata.push_back(xs_f);
298	<	else
299	<	auxdata.push_back(xs_r);
300	<	}
301	<	else {
302	<	if (bh.antisense_align())
303	<	auxdata.push_back(xs_r);
304	<	else
305	<	auxdata.push_back(xs_f);
306	<	}
307	<	}
308	<	else if (library_type == FR_SECONDSTRAND) {
309	<	if (insert_side == FRAG_LEFT \|\| insert_side == FRAG_UNPAIRED){
310	<	if (bh.antisense_align())
311	<	auxdata.push_back(xs_r);
312	<	else
313	<	auxdata.push_back(xs_f);
314	<	}
315	<	else
316	<	{
317	<	if (bh.antisense_align())
318	<	auxdata.push_back(xs_f);
319	<	else
320	<	auxdata.push_back(xs_r);
321	<	}
322	<	}
323	<	} //bh.contiguous()
243	<	if (hitIndex >= 0)
293	>	if (library_type == FR_FIRSTSTRAND) {
294	>	if (insert_side == FRAG_LEFT \|\| insert_side == FRAG_UNPAIRED) {
295	>	if (bh.antisense_align())
296	>	auxdata.push_back(xs_f);
297	>	else
298	>	auxdata.push_back(xs_r);
299	>	}
300	>	else {
301	>	if (bh.antisense_align())
302	>	auxdata.push_back(xs_r);
303	>	else
304	>	auxdata.push_back(xs_f);
305	>	}
306	>	}
307	>	else if (library_type == FR_SECONDSTRAND) {
308	>	if (insert_side == FRAG_LEFT \|\| insert_side == FRAG_UNPAIRED){
309	>	if (bh.antisense_align())
310	>	auxdata.push_back(xs_r);
311	>	else
312	>	auxdata.push_back(xs_f);
313	>	}
314	>	else
315	>	{
316	>	if (bh.antisense_align())
317	>	auxdata.push_back(xs_f);
318	>	else
319	>	auxdata.push_back(xs_r);
320	>	}
321	>	}
322	>	}
323	>	if (hitIndex >= 0)
324		{
325	<	string aux("HI:i:");
326	<	str_appendInt(aux, hitIndex);
327	<	auxdata.push_back(aux);
325	>	string aux("HI:i:");
326	>	str_appendInt(aux, hitIndex);
327	>	auxdata.push_back(aux);
328		}
329		}
330
#	Line 519 \| Line 599
599
600		void print_sam_for_single(const RefSequenceTable& rt,
601		const HitsForRead& hits,
602	<	const FragmentAlignmentGrade& grade,
603	<	FragmentType frag_type,
604	<	//FILE* reads_file,
605	<	Read& read,
606	<	GBamWriter& bam_writer,
607	<	FILE* um_out//write unmapped reads here
528	<	)
602	>	const FragmentAlignmentGrade& grade,
603	>	FragmentType frag_type,
604	>	Read& read,
605	>	GBamWriter& bam_writer,
606	>	FILE* um_out //write unmapped reads here
607	>	)
608		{
609		assert(!read.alt_name.empty());
610	+	if (hits.hits.empty())
611	+	return;
612	+
613		lex_hit_sort s(rt, hits);
614		vector<uint32_t> index_vector;
615
#	Line 537 \| Line 619
619
620		sort(index_vector.begin(), index_vector.end(), s);
621
622	<	size_t primaryHit = (hits.hits.empty()? 0: random() % hits.hits.size());
622	>	size_t primaryHit = 0;
623	>	if (!report_secondary_alignments)
624	>	primaryHit = random() % hits.hits.size();
625	>
626		bool multipleHits = (hits.hits.size() > 1);
627		for (i = 0; i < hits.hits.size(); ++i)
628		{
544	–	bool primary = (i == primaryHit);
629		size_t index = index_vector[i];
630	+	bool primary = (index == primaryHit);
631		const BowtieHit& bh = hits.hits[index];
632		rewrite_sam_record(bam_writer, rt,
633		bh,
#	Line 558 \| Line 643
643		}
644
645		void print_sam_for_pair(const RefSequenceTable& rt,
646	<	const HitsForRead& left_hits,
562	<	const HitsForRead& right_hits,
646	>	const vector<pair<BowtieHit, BowtieHit> >& best_hits,
647		const InsertAlignmentGrade& grade,
648		ReadStream& left_reads_file,
649		ReadStream& right_reads_file,
#	Line 569 \| Line 653
653		uint64_t begin_id = 0,
654		uint64_t end_id = std::numeric_limits<uint64_t>::max())
655		{
572	–	assert (left_hits.insert_id == right_hits.insert_id);
573	–
656		Read left_read;
657		Read right_read;
658	<	assert (left_hits.hits.size() == right_hits.hits.size() \|\|
659	<	(left_hits.hits.empty() \|\| right_hits.hits.empty()));
658	>	if (best_hits.empty())
659	>	return;
660	>
661	>	size_t i;
662	>	HitsForRead right_hits;
663	>	for (i = 0; i < best_hits.size(); ++i)
664	>	right_hits.hits.push_back(best_hits[i].second);
665
666		size_t primaryHit = 0;
667		vector<uint32_t> index_vector;
668	<	if(right_hits.hits.size() > 0)
669	<	{
670	<	lex_hit_sort s(rt, right_hits);
584	<	for (size_t i = 0; i < right_hits.hits.size(); ++i)
585	<	index_vector.push_back(i);
586	<
587	<	sort(index_vector.begin(), index_vector.end(), s);
588	<	primaryHit = random() % right_hits.hits.size();
589	<	}
590	<	else if (left_hits.hits.size() > 0)
591	<	{
592	<	lex_hit_sort s(rt, left_hits);
593	<	for (size_t i = 0; i < left_hits.hits.size(); ++i)
594	<	index_vector.push_back(i);
595	<
596	<	sort(index_vector.begin(), index_vector.end(), s);
597	<	primaryHit = random() % left_hits.hits.size();
598	<	}
668	>	lex_hit_sort s(rt, right_hits);
669	>	for (i = 0; i < right_hits.hits.size(); ++i)
670	>	index_vector.push_back(i);
671
672	<	bool got_left_read = left_reads_file.getRead(left_hits.insert_id, left_read,
601	<	reads_format, false, left_um_out,
602	<	false, begin_id, end_id);
603	<
604	<	bool got_right_read = right_reads_file.getRead(right_hits.insert_id, right_read,
605	<	reads_format, false, right_um_out,
606	<	false, begin_id, end_id);
607	<
608	<	if (left_hits.hits.size() == right_hits.hits.size())
609	<	{
610	<	assert (got_left_read && got_right_read);
611	<	bool multipleHits = (left_hits.hits.size() > 1);
612	<	for (size_t i = 0; i < right_hits.hits.size(); ++i)
613	<	{
614	<	bool primary = (i == primaryHit);
615	<	size_t index = index_vector[i];
616	<	const BowtieHit& right_bh = right_hits.hits[index];
617	<	const BowtieHit& left_bh = left_hits.hits[index];
618	<
619	<	rewrite_sam_record(bam_writer, rt,
620	<	right_bh,
621	<	right_bh.hitfile_rec().c_str(),
622	<	right_read.alt_name.c_str(),
623	<	grade,
624	<	FRAG_RIGHT,
625	<	&left_bh,
626	<	right_hits.hits.size(),
627	<	(i < right_hits.hits.size() - 1) ? &(right_hits.hits[index_vector[i+1]]) : NULL,
628	<	primary,
629	<	(multipleHits? i: -1));
630	<	rewrite_sam_record(bam_writer, rt,
631	<	left_bh,
632	<	left_bh.hitfile_rec().c_str(),
633	<	left_read.alt_name.c_str(),
634	<	grade,
635	<	FRAG_LEFT,
636	<	&right_bh,
637	<	left_hits.hits.size(),
638	<	(i < left_hits.hits.size() - 1) ? &(left_hits.hits[index_vector[i+1]]) : NULL,
639	<	primary,
640	<	(multipleHits? i: -1));
641	<	}
642	<	}
643	<	else if (left_hits.hits.empty())
644	<	{ //only right read was mapped properly
645	<	if (right_um_out) {
646	<	//write it in the mapped file with the #MAPPED# flag
647	<	fprintf(right_um_out, "@%s #MAPPED#\n%s\n+\n%s\n", right_read.alt_name.c_str(),
648	<	right_read.seq.c_str(), right_read.qual.c_str());
649	<	}
650	<	for (size_t i = 0; i < right_hits.hits.size(); ++i)
651	<	{
652	<	bool primary = (i == primaryHit);
653	<	size_t index = index_vector[i];
654	<	const BowtieHit& bh = right_hits.hits[index];
655	<
656	<	rewrite_sam_record(bam_writer, rt,
657	<	bh,
658	<	bh.hitfile_rec().c_str(),
659	<	right_read.alt_name.c_str(),
660	<	grade,
661	<	FRAG_RIGHT,
662	<	NULL,
663	<	right_hits.hits.size(),
664	<	(i < right_hits.hits.size() - 1) ? &(right_hits.hits[index_vector[i+1]]) : NULL,
665	<	primary,
666	<	-1);
667	<
668	<	}
669	<	}
670	<	else if (right_hits.hits.empty())
671	<	{ //only left read was mapped properly
672	<	if (left_um_out) {
673	<	//write it in the mapped file with the #MAPPED# flag
674	<	fprintf(left_um_out, "@%s #MAPPED#\n%s\n+\n%s\n", left_read.alt_name.c_str(), left_read.seq.c_str(),
675	<	left_read.qual.c_str());
676	<	}
672	>	sort(index_vector.begin(), index_vector.end(), s);
673
674	<	for (size_t i = 0; i < left_hits.hits.size(); ++i)
675	<	{
676	<	bool primary = (i == primaryHit);
677	<	size_t index = index_vector[i];
678	<	const BowtieHit& bh = left_hits.hits[index];
679	<	rewrite_sam_record(bam_writer, rt,
680	<	bh,
681	<	bh.hitfile_rec().c_str(),
682	<	left_read.alt_name.c_str(),
683	<	grade,
684	<	FRAG_LEFT,
685	<	NULL,
686	<	left_hits.hits.size(),
687	<	(i < left_hits.hits.size() - 1) ? &(left_hits.hits[index_vector[i+1]]) : NULL,
688	<	primary,
689	<	-1);
690	<
691	<	}
692	<	}
693	<	else
694	<	{
695	<	assert (false);
696	<	}
674	>	if (!report_secondary_alignments)
675	>	primaryHit = random() % right_hits.hits.size();
676	>
677	>	bool got_left_read = left_reads_file.getRead(best_hits[0].first.insert_id(), left_read,
678	>	reads_format, false, begin_id, end_id,
679	>	left_um_out, false);
680	>
681	>	bool got_right_read = right_reads_file.getRead(best_hits[0].first.insert_id(), right_read,
682	>	reads_format, false, begin_id, end_id,
683	>	right_um_out, false);
684	>
685	>	assert (got_left_read && got_right_read);
686	>	bool multipleHits = (best_hits.size() > 1);
687	>	for (i = 0; i < best_hits.size(); ++i)
688	>	{
689	>	size_t index = index_vector[i];
690	>	bool primary = (index == primaryHit);
691	>	const BowtieHit& right_bh = best_hits[index].second;
692	>	const BowtieHit& left_bh = best_hits[index].first;
693	>
694	>	rewrite_sam_record(bam_writer, rt,
695	>	right_bh,
696	>	right_bh.hitfile_rec().c_str(),
697	>	right_read.alt_name.c_str(),
698	>	grade,
699	>	FRAG_RIGHT,
700	>	&left_bh,
701	>	best_hits.size(),
702	>	(i < best_hits.size() - 1) ? &(best_hits[index_vector[i+1]].second) : NULL,
703	>	primary,
704	>	(multipleHits? i: -1));
705	>	rewrite_sam_record(bam_writer, rt,
706	>	left_bh,
707	>	left_bh.hitfile_rec().c_str(),
708	>	left_read.alt_name.c_str(),
709	>	grade,
710	>	FRAG_LEFT,
711	>	&right_bh,
712	>	best_hits.size(),
713	>	(i < best_hits.size() - 1) ? &(best_hits[index_vector[i+1]].first) : NULL,
714	>	primary,
715	>	(multipleHits? i: -1));
716	>	}
717		}
718
719		/**
#	Line 785 \| Line 801
801		for (size_t i = 0; i < hits.hits.size(); ++i)
802		{
803		BowtieHit& bh = hits.hits[i];
804	+	if (bh.edit_dist() > max_read_mismatches)
805	+	continue;
806	+
807		bool filter_hit = false;
808		if (!bh.contiguous())
809		{
#	Line 846 \| Line 865
865		{
866		HitsForRead best_hits;
867		best_hits.insert_id = curr_left_obs_order;
868	<	FragmentAlignmentGrade grade;
868	>	FragmentAlignmentGrade best_grade;
869
870		// Process hits for left singleton, select best alignments
871	<	read_best_alignments(curr_left_hit_group, grade, best_hits, *gtf_junctions);
871	>	read_best_alignments(curr_left_hit_group, best_grade, best_hits, *gtf_junctions);
872		update_junctions(best_hits, *junctions);
873		update_fusions(best_hits, rt, fusions);
874
#	Line 865 \| Line 884
884		{
885		HitsForRead best_hits;
886		best_hits.insert_id = curr_right_obs_order;
887	+	FragmentAlignmentGrade best_grade;
888		if (curr_right_obs_order >= begin_id)
889		{
870	–	FragmentAlignmentGrade grade;
890		// Process hit for right singleton, select best alignments
891	<	read_best_alignments(curr_right_hit_group,grade, best_hits, *gtf_junctions);
891	>	read_best_alignments(curr_right_hit_group, best_grade, best_hits, *gtf_junctions);
892		update_junctions(best_hits, *junctions);
893		update_fusions(best_hits, rt, fusions);
894		}
#	Line 885 \| Line 904
904		curr_left_obs_order < end_id &&
905		curr_left_obs_order != VMAXINT32)
906		{
907	<	if (curr_left_hit_group.hits.empty())
908	<	{
909	<	HitsForRead right_best_hits;
910	<	right_best_hits.insert_id = curr_right_obs_order;
911	<
912	<	FragmentAlignmentGrade grade;
913	<	read_best_alignments(curr_right_hit_group, grade, right_best_hits, *gtf_junctions);
914	<	update_junctions(right_best_hits, *junctions);
915	<	update_fusions(right_best_hits, rt, fusions);
916	<	}
917	<	else if (curr_right_hit_group.hits.empty())
907	>	vector<pair<BowtieHit, BowtieHit> > best_hits;
908	>
909	>	InsertAlignmentGrade grade;
910	>	pair_best_alignments(curr_left_hit_group,
911	>	curr_right_hit_group,
912	>	grade,
913	>	best_hits,
914	>	*gtf_junctions);
915	>
916	>	HitsForRead single_best_hits;
917	>	single_best_hits.insert_id = curr_left_obs_order;
918	>
919	>	// this is too much memory copy - we will make it efficient soon
920	>	for (size_t i = 0; i < best_hits.size(); ++i)
921		{
922	<	HitsForRead left_best_hits;
923	<	left_best_hits.insert_id = curr_left_obs_order;
902	<
903	<	FragmentAlignmentGrade grade;
904	<	// Process hits for left singleton, select best alignments
905	<	read_best_alignments(curr_left_hit_group, grade, left_best_hits, *gtf_junctions);
906	<	update_junctions(left_best_hits, *junctions);
907	<	update_fusions(left_best_hits, rt, fusions);
922	>	single_best_hits.hits.push_back(best_hits[i].first);
923	>	single_best_hits.hits.push_back(best_hits[i].second);
924		}
909	–	else
910	–	{
911	–	HitsForRead left_best_hits;
912	–	HitsForRead right_best_hits;
913	–	left_best_hits.insert_id = curr_left_obs_order;
914	–	right_best_hits.insert_id = curr_right_obs_order;
915	–
916	–	// daehwan - apply gtf_junctions here, too!
925
926	<	InsertAlignmentGrade grade;
927	<	pair_best_alignments(curr_left_hit_group,
928	<	curr_right_hit_group,
921	<	grade,
922	<	left_best_hits,
923	<	right_best_hits);
924	<	update_junctions(left_best_hits, *junctions);
925	<	update_junctions(right_best_hits, *junctions);
926	<	update_fusions(left_best_hits, rt, fusions);
927	<	update_fusions(right_best_hits, rt, fusions);
928	<	}
929	<
926	>	update_junctions(single_best_hits, *junctions);
927	>	update_fusions(single_best_hits, rt, fusions);
928	>
929		l_hs.next_read_hits(curr_left_hit_group);
930		curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id);
931
#	Line 1021 \| Line 1020
1020		best_hits.insert_id = curr_left_obs_order;
1021		FragmentAlignmentGrade grade;
1022		bool got_read=left_reads_file.getRead(curr_left_obs_order, l_read,
1023	<	reads_format, false, left_um_out.file,
1024	<	false, begin_id, end_id);
1025	<	assert(got_read);
1023	>	reads_format, false, begin_id, end_id,
1024	>	left_um_out.file, false);
1025	>	//assert(got_read);
1026
1027		if (right_reads_file.file()) {
1028	+	/*
1029		fprintf(left_um_out.file, "@%s #MAPPED#\n%s\n+\n%s\n", l_read.alt_name.c_str(),
1030		l_read.seq.c_str(), l_read.qual.c_str());
1031	+	*/
1032		got_read=right_reads_file.getRead(curr_left_obs_order, r_read,
1033	<	reads_format, false, right_um_out.file,
1034	<	true, begin_id, end_id);
1035	<	assert(got_read);
1033	>	reads_format, false, begin_id, end_id,
1034	>	right_um_out.file, true);
1035	>	//assert(got_read);
1036		}
1037	<
1037	>
1038		exclude_hits_on_filtered_junctions(*junctions, curr_left_hit_group);
1039
1040		// Process hits for left singleton, select best alignments
1041		read_best_alignments(curr_left_hit_group, grade, best_hits, *gtf_junctions);
1042
1043	<	if (best_hits.hits.size()>0 && best_hits.hits.size() <= max_multihits)
1043	>	if (best_hits.hits.size() > max_multihits)
1044	>	{
1045	>	best_hits.hits.erase(best_hits.hits.begin() + max_multihits, best_hits.hits.end());
1046	>	grade.num_alignments = best_hits.hits.size();
1047	>	}
1048	>
1049	>	if (best_hits.hits.size() > 0)
1050		{
1051		update_junctions(best_hits, *final_junctions);
1052		update_insertions_and_deletions(best_hits, final_insertions, final_deletions);
#	Line 1070 \| Line 1077
1077		if (curr_right_obs_order >= begin_id)
1078		{
1079		bool got_read=right_reads_file.getRead(curr_right_obs_order, r_read,
1080	<	reads_format, false, right_um_out.file,
1081	<	false, begin_id, end_id);
1075	<
1076	<	assert(got_read);
1080	>	reads_format, false, begin_id, end_id,
1081	>	right_um_out.file, false);
1082
1083	+	//assert(got_read);
1084	+	/*
1085		fprintf(right_um_out.file, "@%s #MAPPED#\n%s\n+\n%s\n", r_read.alt_name.c_str(),
1086		r_read.seq.c_str(), r_read.qual.c_str());
1087	+	*/
1088		got_read=left_reads_file.getRead(curr_right_obs_order, l_read,
1089	<	reads_format, false, left_um_out.file,
1090	<	true, begin_id, end_id);
1091	<	assert(got_read);
1089	>	reads_format, false, begin_id, end_id,
1090	>	left_um_out.file, true);
1091	>	//assert(got_read);
1092
1093		exclude_hits_on_filtered_junctions(*junctions, curr_right_hit_group);
1094
1095		// Process hit for right singleton, select best alignments
1096		read_best_alignments(curr_right_hit_group, grade, best_hits, *gtf_junctions);
1097	<	if (best_hits.hits.size()>0 && best_hits.hits.size() <= max_multihits)
1097	>
1098	>	if (best_hits.hits.size() > max_multihits)
1099	>	{
1100	>	best_hits.hits.erase(best_hits.hits.begin() + max_multihits, best_hits.hits.end());
1101	>	grade.num_alignments = best_hits.hits.size();
1102	>	}
1103	>
1104	>	if (best_hits.hits.size() > 0)
1105		{
1106		update_junctions(best_hits, *final_junctions);
1107		update_insertions_and_deletions(best_hits, final_insertions, final_deletions);
#	Line 1120 \| Line 1135
1135		//write it in the mapped file with the #MAPPED# flag
1136
1137		bool got_read=right_reads_file.getRead(curr_left_obs_order, r_read,
1138	<	reads_format, false, right_um_out.file,
1139	<	false, begin_id, end_id);
1140	<	assert(got_read);
1138	>	reads_format, false, begin_id, end_id,
1139	>	right_um_out.file, false);
1140	>	//assert(got_read);
1141	>	/*
1142		fprintf(right_um_out.file, "@%s #MAPPED#\n%s\n+\n%s\n", r_read.alt_name.c_str(),
1143		r_read.seq.c_str(), r_read.qual.c_str());
1144	+	*/
1145		HitsForRead right_best_hits;
1146		right_best_hits.insert_id = curr_right_obs_order;
1147
1148		FragmentAlignmentGrade grade;
1149		read_best_alignments(curr_right_hit_group, grade, right_best_hits, *gtf_junctions);
1150	+
1151	+	if (right_best_hits.hits.size() > max_multihits)
1152	+	{
1153	+	right_best_hits.hits.erase(right_best_hits.hits.begin() + max_multihits, right_best_hits.hits.end());
1154	+	grade.num_alignments = right_best_hits.hits.size();
1155	+	}
1156
1157	<	if (right_best_hits.hits.size()>0 && right_best_hits.hits.size() <= max_multihits)
1157	>	if (right_best_hits.hits.size() > 0)
1158		{
1159		update_junctions(right_best_hits, *final_junctions);
1160		update_insertions_and_deletions(right_best_hits, final_insertions, final_deletions);
#	Line 1152 \| Line 1175
1175		left_best_hits.insert_id = curr_left_obs_order;
1176		//only left read mapped
1177		bool got_read=left_reads_file.getRead(curr_left_obs_order, l_read,
1178	<	reads_format, false, left_um_out.file,
1179	<	false, begin_id, end_id);
1180	<	assert(got_read);
1178	>	reads_format, false, begin_id, end_id,
1179	>	left_um_out.file, false);
1180	>	//assert(got_read);
1181	>	/*
1182		fprintf(left_um_out.file, "@%s #MAPPED#\n%s\n+\n%s\n", l_read.alt_name.c_str(),
1183		l_read.seq.c_str(), l_read.qual.c_str());
1184	+	*/
1185		FragmentAlignmentGrade grade;
1186		// Process hits for left singleton, select best alignments
1187		read_best_alignments(curr_left_hit_group, grade, left_best_hits, *gtf_junctions);
1188	<	if (left_best_hits.hits.size()>0 && left_best_hits.hits.size() <= max_multihits)
1188	>
1189	>	if (left_best_hits.hits.size() > max_multihits)
1190	>	{
1191	>	left_best_hits.hits.erase(left_best_hits.hits.begin() + max_multihits, left_best_hits.hits.end());
1192	>	grade.num_alignments = left_best_hits.hits.size();
1193	>	}
1194	>
1195	>	if (left_best_hits.hits.size() > 0)
1196		{
1197		update_junctions(left_best_hits, *final_junctions);
1198		update_insertions_and_deletions(left_best_hits, final_insertions, final_deletions);
#	Line 1176 \| Line 1208
1208		}
1209		}
1210		else
1211	<	{ //hits for both left and right reads
1212	<	HitsForRead left_best_hits;
1181	<	HitsForRead right_best_hits;
1182	<	left_best_hits.insert_id = curr_left_obs_order;
1183	<	right_best_hits.insert_id = curr_right_obs_order;
1211	>	{
1212	>	vector<pair<BowtieHit, BowtieHit> > best_hits;
1213
1214		InsertAlignmentGrade grade;
1215		pair_best_alignments(curr_left_hit_group,
1216		curr_right_hit_group,
1217		grade,
1218	<	left_best_hits,
1219	<	right_best_hits);
1218	>	best_hits,
1219	>	*gtf_junctions);
1220	>
1221	>	if (best_hits.size() > max_multihits)
1222	>	{
1223	>	best_hits.erase(best_hits.begin() + max_multihits, best_hits.end());
1224	>	grade.num_alignments = best_hits.size();
1225	>	}
1226	>
1227	>	//hits for both left and right reads
1228	>	HitsForRead single_best_hits;
1229	>	single_best_hits.insert_id = curr_left_obs_order;
1230	>	for (size_t i = 0; i < best_hits.size(); ++i)
1231	>	{
1232	>	single_best_hits.hits.push_back(best_hits[i].first);
1233	>	single_best_hits.hits.push_back(best_hits[i].second);
1234	>	}
1235
1236	<	if (left_best_hits.hits.size()>0 && left_best_hits.hits.size() <= max_multihits &&
1193	<	right_best_hits.hits.size()>0 && right_best_hits.hits.size() <= max_multihits)
1236	>	if (single_best_hits.hits.size() > 0)
1237		{
1238	<	update_junctions(left_best_hits, *final_junctions);
1239	<	update_junctions(right_best_hits, *final_junctions);
1197	<	update_insertions_and_deletions(left_best_hits, final_insertions, final_deletions);
1198	<	update_insertions_and_deletions(right_best_hits, final_insertions, final_deletions);
1238	>	update_junctions(single_best_hits, *final_junctions);
1239	>	update_insertions_and_deletions(single_best_hits, final_insertions, final_deletions);
1240
1241	<	pair_support(left_best_hits, right_best_hits, final_fusions, fusions);
1242	<	update_fusions(left_best_hits, rt, final_fusions, *fusions);
1202	<	update_fusions(right_best_hits, rt, final_fusions, *fusions);
1241	>	pair_support(best_hits, final_fusions, fusions);
1242	>	update_fusions(single_best_hits, rt, final_fusions, *fusions);
1243
1244		print_sam_for_pair(*rt,
1245	<	left_best_hits,
1206	<	right_best_hits,
1245	>	best_hits,
1246		grade,
1247		left_reads_file,
1248		right_reads_file,
#	Line 1224 \| Line 1263
1263		} //while we still have unreported hits..
1264		//print the remaining unmapped reads at the end of each reads' stream
1265		left_reads_file.getRead(VMAXINT32, l_read,
1266	<	reads_format, false, left_um_out.file,
1267	<	false, begin_id, end_id);
1266	>	reads_format, false, begin_id, end_id,
1267	>	left_um_out.file, false);
1268		if (right_reads_file.file())
1269		right_reads_file.getRead(VMAXINT32, r_read,
1270	<	reads_format, false, right_um_out.file,
1271	<	false, begin_id, end_id);
1270	>	reads_format, false, begin_id, end_id,
1271	>	right_um_out.file, false);
1272
1273		left_um_out.close();
1274		right_um_out.close();

Diff Legend

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