[ViewVC] Diff of: gclib/tophat_cpp/tophat

#	Line 31 \| Line 31
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"
#	Line 74 \| Line 75
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		}
#	Line 144 \| Line 147
147		}
148		}
149
150	<	bool set_insert_alignment_grade(const BowtieHit& lh, const BowtieHit& rh, InsertAlignmentGrade& grade)
150	>	bool set_insert_alignment_grade(const BowtieHit& lh, const BowtieHit& rh, const JunctionSet& junctions, InsertAlignmentGrade& grade)
151		{
149	–	// max mate inner distance (genomic)
150	–	int min_mate_inner_dist = inner_dist_mean - inner_dist_std_dev;
151	–	if (max_mate_inner_dist == -1)
152	–	max_mate_inner_dist = inner_dist_mean + inner_dist_std_dev;
153	–
152		bool fusion = false;
153		bool left_fusion = lh.fusion_opcode() != FUSION_NOTHING;
154		bool right_fusion = rh.fusion_opcode() != FUSION_NOTHING;
#	Line 180 \| Line 178
178		}
179		}
180
181	<	grade = InsertAlignmentGrade(lh, rh, min_mate_inner_dist, max_mate_inner_dist, fusion);
181	>	grade = InsertAlignmentGrade(lh, rh, junctions, fusion);
182
183		return true;
184		}
185
186		struct cmp_pair_alignment
187		{
188	<	cmp_pair_alignment(const JunctionSet& gtf_juncs) :
189	<	_gtf_juncs(&gtf_juncs) {}
188	>	cmp_pair_alignment(const JunctionSet& junctions) :
189	>	_junctions(&junctions) {}
190
191		bool operator() (const pair<BowtieHit, BowtieHit>& l, const pair<BowtieHit, BowtieHit>& r) const
192		{
193	<	InsertAlignmentGrade gl; set_insert_alignment_grade(l.first, l.second, gl);
194	<	InsertAlignmentGrade gr; set_insert_alignment_grade(r.first, r.second, gr);
193	>	InsertAlignmentGrade gl; set_insert_alignment_grade(l.first, l.second, *_junctions, gl);
194	>	InsertAlignmentGrade gr; set_insert_alignment_grade(r.first, r.second, *_junctions, gr);
195
196		bool better = gr < gl;
197		bool worse = gl < gr;
#	Line 204 \| Line 202
202		return false;
203		}
204
205	<	const JunctionSet* _gtf_juncs;
205	>	const JunctionSet* _junctions;
206		};
207
208		void pair_best_alignments(const HitsForRead& left_hits,
#	Line 243 \| Line 241
241		junctions, insertions, deletions, fusions, coverage);
242		rh.alignment_score(align_status._alignment_score);
243		InsertAlignmentGrade g;
244	<	bool allowed = set_insert_alignment_grade(lh, rh, g);
244	>	bool allowed;
245	>	allowed = set_insert_alignment_grade(lh, rh, final_report ? junctions : gtf_junctions, g);
246	>
247	>	// daehwan - for debugging purposes
248	>	if (lh.insert_id() == 325708 && !g.is_fusion() && false)
249	>	{
250	>	fprintf(stderr, "lh %d:%d %s score: %d (from %d) NM: %d\n",
251	>	lh.ref_id(), lh.left(), print_cigar(lh.cigar()).c_str(),
252	>	lh.alignment_score(), left[i].alignment_score(), lh.edit_dist());
253	>	fprintf(stderr, "rh %d:%d %s score: %d (from %d) NM: %d\n",
254	>	rh.ref_id(), rh.left(), print_cigar(rh.cigar()).c_str(),
255	>	rh.alignment_score(), right[j].alignment_score(), rh.edit_dist());
256	>	fprintf(stderr, "combined score: %d is_fusion(%d)\n", g.align_score(), g.is_fusion());
257	>	}
258
259		if (!allowed) continue;
260		if (!fusion_search && !report_discordant_pair_alignments && g.is_fusion()) continue;
#	Line 271 \| Line 282
282
283		if ((report_secondary_alignments \|\| !final_report) && best_hits.size() > 0)
284		{
285	<	cmp_pair_alignment cmp(gtf_junctions);
285	>	cmp_pair_alignment cmp(final_report ? junctions : gtf_junctions);
286		sort(best_hits.begin(), best_hits.end(), cmp);
287	<	set_insert_alignment_grade(best_hits[0].first, best_hits[0].second, best_grade);
287	>	set_insert_alignment_grade(best_hits[0].first, best_hits[0].second, final_report ? junctions : gtf_junctions, best_grade);
288		}
289
290		if (final_report)
#	Line 382 \| Line 393
393		string ref_name = sam_toks[2], ref_name2 = "";
394		char cigar1[1024] = {0}, cigar2[1024] = {0};
395		string left_seq, right_seq, left_qual, right_qual;
396	<	int left = -1, left1 = -1, left2 = -1;
396	>	int left1 = -1, left2 = -1;
397		bool fusion_alignment = false;
398		size_t XF_index = 0;
399		for (size_t i = 11; i < sam_toks.size(); ++i)
#	Line 413 \| Line 424
424		else if (strncmp(tok.c_str(), "AS", 2) == 0)
425		{
426		char AS_score[128] = {0};
427	<	sprintf(AS_score, "AS:i:%d", bh.alignment_score());
427	>	sprintf(AS_score, "AS:i:%d", min(0, bh.alignment_score()));
428		tok = AS_score;
429		}
430		}
#	Line 441 \| Line 452
452		flag \|= 0x100;
453
454		int gpos=isdigit(sam_toks[3][0]) ? atoi(sam_toks[3].c_str()) : 0;
455	<	int mapQ=255;
455	>	int mapQ = 50;
456		if (num_hits > 1) {
457		double err_prob = 1 - (1.0 / num_hits);
458		mapQ = (int)(-10.0 * log(err_prob) / log(10.0));
#	Line 449 \| Line 460
460		int tlen =atoi(sam_toks[8].c_str()); //TLEN
461		int mate_pos=atoi(sam_toks[7].c_str());
462
463	+	string rg_aux = "";
464	+	if (!sam_readgroup_id.empty())
465	+	rg_aux = string("RG:Z:") + sam_readgroup_id;
466	+
467		GBamRecord* bamrec=NULL;
468		if (fusion_alignment) {
469		vector<string> auxdata;
470		add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex);
471	+	if (rg_aux != "")
472	+	auxdata.push_back(rg_aux);
473		bamrec=bam_writer.new_record(qname.c_str(), flag, ref_name.c_str(), left1 + 1, mapQ,
474		cigar1, sam_toks[6].c_str(), mate_pos,
475		tlen, left_seq.c_str(), left_qual.c_str(), &auxdata);
#	Line 462 \| Line 479
479		auxdata.clear();
480		sam_toks[XF_index][5] = '2';
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_name2.c_str(), left2 + 1, mapQ,
485		cigar2, sam_toks[6].c_str(), mate_pos,
486		tlen, right_seq.c_str(), right_qual.c_str(), &auxdata);
#	Line 470 \| Line 489
489		} else {
490		vector<string> auxdata;
491		add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex);
492	+	if (rg_aux != "")
493	+	auxdata.push_back(rg_aux);
494		bamrec=bam_writer.new_record(qname.c_str(), flag, sam_toks[2].c_str(), gpos, mapQ,
495		sam_toks[5].c_str(), sam_toks[6].c_str(), mate_pos,
496		tlen, sam_toks[9].c_str(), sam_toks[10].c_str(), &auxdata);
#	Line 516 \| Line 537
537		string ref_name = sam_toks[2], ref_name2 = "";
538		char cigar1[1024] = {0}, cigar2[1024] = {0};
539		string left_seq, right_seq, left_qual, right_qual;
540	<	int left = -1, left1 = -1, left2 = -1;
540	>	int left1 = -1, left2 = -1;
541		bool fusion_alignment = false;
542		size_t XF_tok_idx = 11;
543		for (; XF_tok_idx < sam_toks.size(); ++XF_tok_idx)
#	Line 546 \| Line 567
567		else if (strncmp(tok.c_str(), "AS", 2) == 0)
568		{
569		char AS_score[128] = {0};
570	<	sprintf(AS_score, "AS:i:%d", bh.alignment_score());
570	>	sprintf(AS_score, "AS:i:%d", min(0, bh.alignment_score()));
571		tok = AS_score;
572		}
573		}
574
575		int gpos=isdigit(sam_toks[3][0]) ? atoi(sam_toks[3].c_str()) : 0;
576	<	int mapQ=255;
576	>	int mapQ = 50;
577		if (grade.num_alignments > 1) {
578		double err_prob = 1 - (1.0 / grade.num_alignments);
579		mapQ = (int)(-10.0 * log(err_prob) / log(10.0));
#	Line 589 \| Line 610
610		flag \|= 0x0008;
611		}
612
613	+	string rg_aux = "";
614	+	if (!sam_readgroup_id.empty())
615	+	rg_aux = string("RG:Z:") + sam_readgroup_id;
616	+
617		GBamRecord* bamrec=NULL;
618		if (fusion_alignment) {
619		vector<string> auxdata;
620		add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex);
621	+	if (rg_aux != "")
622	+	auxdata.push_back(rg_aux);
623		bamrec=bam_writer.new_record(qname.c_str(), flag, ref_name.c_str(), left1 + 1, mapQ,
624		cigar1, sam_toks[6].c_str(), mate_pos,
625		tlen, left_seq.c_str(), left_qual.c_str(), &auxdata);
#	Line 602 \| Line 629
629		auxdata.clear();
630		sam_toks[XF_tok_idx][5] = '2';
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_name2.c_str(), left2 + 1, mapQ,
635		cigar2, sam_toks[6].c_str(), mate_pos,
636		tlen, right_seq.c_str(), right_qual.c_str(), &auxdata);
#	Line 610 \| Line 639
639		} else {
640		vector<string> auxdata;
641		add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex);
642	+	if (rg_aux != "")
643	+	auxdata.push_back(rg_aux);
644		bamrec=bam_writer.new_record(qname.c_str(), flag, sam_toks[2].c_str(), gpos, mapQ,
645		sam_toks[5].c_str(), sam_toks[6].c_str(), mate_pos,
646		tlen, sam_toks[9].c_str(), sam_toks[10].c_str(), &auxdata);
#	Line 646 \| Line 677
677		void print_sam_for_single(const RefSequenceTable& rt,
678		const HitsForRead& hits,
679		FragmentType frag_type,
680	<	Read& read,
680	>	const string& alt_name,
681		GBamWriter& bam_writer,
682	<	GBamWriter* um_out //write unmapped reads here
652	<	)
682	>	boost::random::mt19937& rng)
683		{
684	<	assert(!read.alt_name.empty());
684	>	//assert(!read.alt_name.empty());
685		if (hits.hits.empty())
686		return;
687
#	Line 666 \| Line 696
696
697		size_t primaryHit = 0;
698		if (!report_secondary_alignments)
699	<	primaryHit = random() % hits.hits.size();
699	>	primaryHit = rng() % hits.hits.size();
700
701		bool multipleHits = (hits.hits.size() > 1);
702		for (i = 0; i < hits.hits.size(); ++i)
#	Line 677 \| Line 707
707		rewrite_sam_record(bam_writer, rt,
708		bh,
709		bh.hitfile_rec().c_str(),
710	<	read.alt_name.c_str(),
710	>	alt_name.c_str(),
711		frag_type,
712		hits.hits.size(),
713		(i < hits.hits.size()-1) ? &(hits.hits[index_vector[i+1]]) : NULL,
#	Line 689 \| Line 719
719		void print_sam_for_pair(const RefSequenceTable& rt,
720		const vector<pair<BowtieHit, BowtieHit> >& best_hits,
721		const InsertAlignmentGrade& grade,
692	–	ReadStream& left_reads_file,
693	–	ReadStream& right_reads_file,
722		GBamWriter& bam_writer,
723	<	GBamWriter* left_um_out,
724	<	GBamWriter* right_um_out,
723	>	const string& left_alt_name,
724	>	const string& right_alt_name,
725	>	boost::random::mt19937& rng,
726		uint64_t begin_id = 0,
727		uint64_t end_id = std::numeric_limits<uint64_t>::max())
728		{
#	Line 716 \| Line 745
745		sort(index_vector.begin(), index_vector.end(), s);
746
747		if (!report_secondary_alignments)
748	<	primaryHit = random() % right_hits.hits.size();
748	>	primaryHit = rng() % right_hits.hits.size();
749
721	–	bool got_left_read = left_reads_file.getRead(best_hits[0].first.insert_id(), left_read,
722	–	reads_format, false, begin_id, end_id,
723	–	left_um_out, false);
724	–
725	–	bool got_right_read = right_reads_file.getRead(best_hits[0].first.insert_id(), right_read,
726	–	reads_format, false, begin_id, end_id,
727	–	right_um_out, false);
728	–
729	–	assert (got_left_read && got_right_read);
750		bool multipleHits = (best_hits.size() > 1);
751		for (i = 0; i < best_hits.size(); ++i)
752		{
#	Line 738 \| Line 758
758		rewrite_sam_record(bam_writer, rt,
759		right_bh,
760		right_bh.hitfile_rec().c_str(),
761	<	right_read.alt_name.c_str(),
761	>	right_alt_name.c_str(),
762		grade,
763		FRAG_RIGHT,
764		&left_bh,
#	Line 749 \| Line 769
769		rewrite_sam_record(bam_writer, rt,
770		left_bh,
771		left_bh.hitfile_rec().c_str(),
772	<	left_read.alt_name.c_str(),
772	>	left_alt_name.c_str(),
773		grade,
774		FRAG_LEFT,
775		&right_bh,
#	Line 898 \| Line 918
918		hits = remaining;
919		}
920
921	+	void realign_reads(HitsForRead& hits,
922	+	const RefSequenceTable& rt,
923	+	const JunctionSet& junctions,
924	+	const JunctionSet& rev_junctions,
925	+	const InsertionSet& insertions,
926	+	const DeletionSet& deletions,
927	+	const DeletionSet& rev_deletions,
928	+	const FusionSet& fusions)
929	+	{
930	+	vector<BowtieHit> additional_hits;
931	+	for (size_t i = 0; i < hits.hits.size(); ++i)
932	+	{
933	+	BowtieHit& bh = hits.hits[i];
934	+
935	+	const vector<CigarOp>& cigars = bh.cigar();
936	+	int pos = bh.left();
937	+	int refid = bh.ref_id();
938	+	for (size_t j = 0; j < cigars.size(); ++j)
939	+	{
940	+	const CigarOp& op = cigars[j];
941	+	if (j == 0 \|\| j == cigars.size() - 1)
942	+	{
943	+	// let's do this for MATCH case only,
944	+	if (op.opcode == MATCH \|\| op.opcode == mATCH)
945	+	{
946	+	int left1, left2;
947	+	if (op.opcode == MATCH)
948	+	{
949	+	left1 = pos;
950	+	left2 = pos + op.length - 1;
951	+	}
952	+	else
953	+	{
954	+	left1 = pos - op.length + 1;
955	+	left2 = pos;
956	+	}
957	+
958	+	{
959	+	JunctionSet::const_iterator lb, ub;
960	+	JunctionSet temp_junctions;
961	+	if (j == 0)
962	+	{
963	+	lb = rev_junctions.upper_bound(Junction(refid, left1, 0, true));
964	+	ub = rev_junctions.lower_bound(Junction(refid, left2, left2, false));
965	+	while (lb != ub && lb != rev_junctions.end())
966	+	{
967	+	Junction temp_junction = lb->first;
968	+	temp_junction.left = lb->first.right;
969	+	temp_junction.right = lb->first.left;
970	+	temp_junctions[temp_junction] = lb->second;
971	+	++lb;
972	+	}
973	+	}
974	+
975	+	if (j == cigars.size() - 1)
976	+	{
977	+	int common_right = left2 + max_report_intron_length;
978	+	lb = junctions.upper_bound(Junction(refid, left1, common_right, true));
979	+	ub = junctions.lower_bound(Junction(refid, left2, common_right, false));
980	+	while (lb != ub && lb != junctions.end())
981	+	{
982	+	temp_junctions[lb->first] = lb->second;
983	+	++lb;
984	+	}
985	+	}
986	+
987	+	JunctionSet::const_iterator junc_iter = temp_junctions.begin();
988	+	for (; junc_iter != temp_junctions.end(); ++junc_iter)
989	+	{
990	+	Junction junc = junc_iter->first;
991	+
992	+	/*
993	+	fprintf(stderr, "(type%d) %d %d-%d %s (AS:%d XM:%d) with junc %u-%u\n",
994	+	!use_rev_junctions,
995	+	bh.insert_id(), bh.left(), bh.right(),
996	+	print_cigar(bh.cigar()).c_str(),
997	+	bh.alignment_score(), bh.edit_dist(),
998	+	junc.left, junc.right);
999	+	//*/
1000	+
1001	+	int new_left = bh.left();
1002	+	int intron_length = junc.right - junc.left - 1;
1003	+	vector<CigarOp> new_cigars;
1004	+	bool anchored = false;
1005	+	if (j == 0 && bh.left() > (int)junc.left)
1006	+	{
1007	+	new_left -= intron_length;
1008	+	int before_match_length = junc.left - new_left + 1;;
1009	+	int after_match_length = op.length - before_match_length;
1010	+
1011	+	if (before_match_length > 0 && after_match_length)
1012	+	{
1013	+	anchored = true;
1014	+	new_cigars.push_back(CigarOp(MATCH, before_match_length));
1015	+	new_cigars.push_back(CigarOp(REF_SKIP, intron_length));
1016	+	new_cigars.push_back(CigarOp(MATCH, after_match_length));
1017	+
1018	+	new_cigars.insert(new_cigars.end(), cigars.begin() + 1, cigars.end());
1019	+	}
1020	+	}
1021	+	else
1022	+	{
1023	+	new_cigars.insert(new_cigars.end(), cigars.begin(), cigars.end() - 1);
1024	+
1025	+	int before_match_length = junc.left - pos + 1;
1026	+	int after_match_length = op.length - before_match_length;
1027	+
1028	+	if (before_match_length > 0 && after_match_length > 0)
1029	+	{
1030	+	anchored = true;
1031	+	new_cigars.push_back(CigarOp(MATCH, before_match_length));
1032	+	new_cigars.push_back(CigarOp(REF_SKIP, intron_length));
1033	+	new_cigars.push_back(CigarOp(MATCH, after_match_length));
1034	+	}
1035	+	}
1036	+
1037	+	if (anchored)
1038	+	{
1039	+	BowtieHit new_bh(bh.ref_id(),
1040	+	bh.ref_id2(),
1041	+	bh.insert_id(),
1042	+	new_left,
1043	+	new_cigars,
1044	+	bh.antisense_align(),
1045	+	junc.antisense,
1046	+	0, /* edit_dist - needs to be recalculated */
1047	+	0, /* splice_mms - needs to be recalculated */
1048	+	false);
1049	+
1050	+	new_bh.seq(bh.seq());
1051	+	new_bh.qual(bh.qual());
1052	+
1053	+	vector<string> aux_fields;
1054	+	bowtie_sam_extra(new_bh, rt, aux_fields);
1055	+
1056	+	vector<string>::const_iterator aux_iter = aux_fields.begin();
1057	+	for (; aux_iter != aux_fields.end(); ++aux_iter)
1058	+	{
1059	+	const string& aux_field = *aux_iter;
1060	+	if (strncmp(aux_field.c_str(), "AS", 2) == 0)
1061	+	{
1062	+	int alignment_score = atoi(aux_field.c_str() + 5);
1063	+	new_bh.alignment_score(alignment_score);
1064	+	}
1065	+	else if (strncmp(aux_field.c_str(), "XM", 2) == 0)
1066	+	{
1067	+	int XM_value = atoi(aux_field.c_str() + 5);
1068	+	new_bh.edit_dist(XM_value);
1069	+	}
1070	+	}
1071	+
1072	+	vector<string> sam_toks;
1073	+	tokenize(bh.hitfile_rec().c_str(), "\t", sam_toks);
1074	+
1075	+	char coord[20] = {0,};
1076	+	sprintf(coord, "%d", new_bh.left() + 1);
1077	+	sam_toks[3] = coord;
1078	+	sam_toks[5] = print_cigar(new_bh.cigar());
1079	+	for (size_t a = 11; a < sam_toks.size(); ++a)
1080	+	{
1081	+	string& sam_tok = sam_toks[a];
1082	+	for (size_t b = 0; b < aux_fields.size(); ++b)
1083	+	{
1084	+	const string& aux_tok = aux_fields[b];
1085	+	if (strncmp(sam_tok.c_str(), aux_tok.c_str(), 5) == 0)
1086	+	{
1087	+	sam_tok = aux_tok;
1088	+	break;
1089	+	}
1090	+	}
1091	+	}
1092	+
1093	+	if (!bh.is_spliced())
1094	+	{
1095	+	if (junc.antisense)
1096	+	sam_toks.push_back("XS:A:-");
1097	+	else
1098	+	sam_toks.push_back("XS:A:+");
1099	+	}
1100	+
1101	+
1102	+	string new_rec = "";
1103	+	for (size_t d = 0; d < sam_toks.size(); ++d)
1104	+	{
1105	+	new_rec += sam_toks[d];
1106	+	if (d < sam_toks.size() - 1)
1107	+	new_rec += "\t";
1108	+	}
1109	+
1110	+	new_bh.hitfile_rec(new_rec);
1111	+
1112	+	if (new_bh.edit_dist() <= bh.edit_dist())
1113	+	additional_hits.push_back(new_bh);
1114	+
1115	+	/*
1116	+	fprintf(stderr, "\t%d %d-%d %s (AS:%d XM:%d) with junc %u-%u\n",
1117	+	new_bh.insert_id(), new_bh.left(), new_bh.right(),
1118	+	print_cigar(new_bh.cigar()).c_str(),
1119	+	new_bh.alignment_score(), new_bh.edit_dist(),
1120	+	junc.left, junc.right);
1121	+	//*/
1122	+	}
1123	+	}
1124	+	}
1125	+
1126	+	#if 0
1127	+	{
1128	+	DeletionSet::const_iterator lb, ub;
1129	+	bool use_rev_deletions = (j == 0);
1130	+	const DeletionSet& curr_deletions = (use_rev_deletions ? rev_deletions : deletions);
1131	+	if (use_rev_deletions)
1132	+	{
1133	+	lb = curr_deletions.upper_bound(Deletion(refid, left1, 0, true));
1134	+	ub = curr_deletions.lower_bound(Deletion(refid, left2, left2, false));
1135	+	}
1136	+	else
1137	+	{
1138	+	int common_right = left2 + 100;
1139	+	lb = curr_deletions.upper_bound(Deletion(refid, left1, common_right, true));
1140	+	ub = curr_deletions.lower_bound(Deletion(refid, left2, common_right, false));
1141	+	}
1142	+
1143	+	while (lb != curr_deletions.end() && lb != ub)
1144	+	{
1145	+	Deletion del = lb->first;
1146	+	if (use_rev_deletions)
1147	+	{
1148	+	int temp = del.left;
1149	+	del.left = del.right;
1150	+	del.right = temp;
1151	+	}
1152	+
1153	+	// daehwan - for debuggin purposes
1154	+	/*
1155	+	fprintf(stderr, "(type%d) %d %d-%d %s (AS:%d XM:%d) with junc %u-%u\n",
1156	+	!use_rev_junctions,
1157	+	bh.insert_id(), bh.left(), bh.right(),
1158	+	print_cigar(bh.cigar()).c_str(),
1159	+	bh.alignment_score(), bh.edit_dist(),
1160	+	junc.left, junc.right);
1161	+	*/
1162	+
1163	+	int del_length = del.right - del.left - 1;
1164	+	int new_left = bh.left();
1165	+	if (j == 0)
1166	+	new_left -= del_length;
1167	+
1168	+	vector<CigarOp> new_cigars;
1169	+	if (j == 0)
1170	+	{
1171	+	int before_match_length = del.left - new_left + 1;;
1172	+	int after_match_length = op.length - before_match_length;
1173	+
1174	+	if (before_match_length > 0)
1175	+	new_cigars.push_back(CigarOp(MATCH, before_match_length));
1176	+	new_cigars.push_back(CigarOp(DEL, del_length));
1177	+	if (after_match_length > 0)
1178	+	new_cigars.push_back(CigarOp(MATCH, after_match_length));
1179	+
1180	+	new_cigars.insert(new_cigars.end(), cigars.begin() + 1, cigars.end());
1181	+	}
1182	+	else
1183	+	{
1184	+	new_cigars.insert(new_cigars.end(), cigars.begin(), cigars.end() - 1);
1185	+
1186	+	int before_match_length = del.left - pos + 1;
1187	+	int after_match_length = op.length - before_match_length;
1188	+
1189	+	if (before_match_length > 0)
1190	+	new_cigars.push_back(CigarOp(MATCH, before_match_length));
1191	+	new_cigars.push_back(CigarOp(DEL, del_length));
1192	+	if (after_match_length > 0)
1193	+	new_cigars.push_back(CigarOp(MATCH, after_match_length));
1194	+	}
1195	+
1196	+	BowtieHit new_bh(bh.ref_id(),
1197	+	bh.ref_id2(),
1198	+	bh.insert_id(),
1199	+	new_left,
1200	+	new_cigars,
1201	+	bh.antisense_align(),
1202	+	bh.antisense_splice(),
1203	+	0, /* edit_dist - needs to be recalculated */
1204	+	0, /* splice_mms - needs to be recalculated */
1205	+	false);
1206	+
1207	+	new_bh.seq(bh.seq());
1208	+	new_bh.qual(bh.qual());
1209	+
1210	+	vector<string> aux_fields;
1211	+	bowtie_sam_extra(new_bh, rt, aux_fields);
1212	+
1213	+	vector<string>::const_iterator aux_iter = aux_fields.begin();
1214	+	for (; aux_iter != aux_fields.end(); ++aux_iter)
1215	+	{
1216	+	const string& aux_field = *aux_iter;
1217	+	if (strncmp(aux_field.c_str(), "AS", 2) == 0)
1218	+	{
1219	+	int alignment_score = atoi(aux_field.c_str() + 5);
1220	+	new_bh.alignment_score(alignment_score);
1221	+	}
1222	+	else if (strncmp(aux_field.c_str(), "XM", 2) == 0)
1223	+	{
1224	+	int XM_value = atoi(aux_field.c_str() + 5);
1225	+	new_bh.edit_dist(XM_value);
1226	+	}
1227	+	}
1228	+
1229	+	vector<string> sam_toks;
1230	+	tokenize(bh.hitfile_rec().c_str(), "\t", sam_toks);
1231	+
1232	+	char coord[20] = {0,};
1233	+	sprintf(coord, "%d", new_bh.left() + 1);
1234	+	sam_toks[3] = coord;
1235	+	sam_toks[5] = print_cigar(new_bh.cigar());
1236	+	for (size_t a = 11; a < sam_toks.size(); ++a)
1237	+	{
1238	+	string& sam_tok = sam_toks[a];
1239	+	for (size_t b = 0; b < aux_fields.size(); ++b)
1240	+	{
1241	+	const string& aux_tok = aux_fields[b];
1242	+	if (strncmp(sam_tok.c_str(), aux_tok.c_str(), 5) == 0)
1243	+	{
1244	+	sam_tok = aux_tok;
1245	+	break;
1246	+	}
1247	+	}
1248	+	}
1249	+
1250	+	string new_rec = "";
1251	+	for (size_t d = 0; d < sam_toks.size(); ++d)
1252	+	{
1253	+	new_rec += sam_toks[d];
1254	+	if (d < sam_toks.size() - 1)
1255	+	new_rec += "\t";
1256	+	}
1257	+
1258	+	new_bh.hitfile_rec(new_rec);
1259	+
1260	+	if (new_bh.edit_dist() <= bh.edit_dist())
1261	+	additional_hits.push_back(new_bh);
1262	+
1263	+	/*
1264	+	fprintf(stderr, "\t%d %d-%d %s (AS:%d XM:%d) with junc %u-%u\n",
1265	+	new_bh.insert_id(), new_bh.left(), new_bh.right(),
1266	+	print_cigar(new_bh.cigar()).c_str(),
1267	+	new_bh.alignment_score(), new_bh.edit_dist(),
1268	+	junc.left, junc.right);
1269	+	*/
1270	+
1271	+	++lb;
1272	+	}
1273	+	}
1274	+
1275	+	{
1276	+	InsertionSet::const_iterator lb, ub;
1277	+	lb = insertions.upper_bound(Insertion(refid, left1, ""));
1278	+	ub = insertions.lower_bound(Insertion(refid, left2, ""));
1279	+
1280	+	while (lb != insertions.end() && lb != ub)
1281	+	{
1282	+	// daehwan - for debugging purposse
1283	+	break;
1284	+
1285	+	Insertion ins = lb->first;
1286	+
1287	+	// daehwan - for debugging purposes
1288	+	/*
1289	+	fprintf(stderr, "(type%d) %d %d-%d %s (AS:%d XM:%d) with junc %u-%u\n",
1290	+	!use_rev_junctions,
1291	+	bh.insert_id(), bh.left(), bh.right(),
1292	+	print_cigar(bh.cigar()).c_str(),
1293	+	bh.alignment_score(), bh.edit_dist(),
1294	+	junc.left, junc.right);
1295	+	*/
1296	+
1297	+	int ins_length = ins.sequence.length();
1298	+	int new_left = bh.left();
1299	+	if (j == 0)
1300	+	new_left -= ins_length;
1301	+
1302	+	vector<CigarOp> new_cigars;
1303	+	if (j == 0)
1304	+	{
1305	+	int before_match_length = ins.left - new_left + 1;;
1306	+	int after_match_length = op.length - before_match_length - ins_length;
1307	+
1308	+	if (before_match_length > 0)
1309	+	new_cigars.push_back(CigarOp(MATCH, before_match_length));
1310	+	new_cigars.push_back(CigarOp(INS, ins_length));
1311	+	if (after_match_length > 0)
1312	+	new_cigars.push_back(CigarOp(MATCH, after_match_length));
1313	+
1314	+	new_cigars.insert(new_cigars.end(), cigars.begin() + 1, cigars.end());
1315	+	}
1316	+	else
1317	+	{
1318	+	new_cigars.insert(new_cigars.end(), cigars.begin(), cigars.end() - 1);
1319	+
1320	+	int before_match_length = ins.left - pos + 1;
1321	+	int after_match_length = op.length - before_match_length - ins_length;
1322	+
1323	+	if (before_match_length > 0)
1324	+	new_cigars.push_back(CigarOp(MATCH, before_match_length));
1325	+	new_cigars.push_back(CigarOp(INS, ins_length));
1326	+	if (after_match_length > 0)
1327	+	new_cigars.push_back(CigarOp(MATCH, after_match_length));
1328	+	}
1329	+
1330	+	BowtieHit new_bh(bh.ref_id(),
1331	+	bh.ref_id2(),
1332	+	bh.insert_id(),
1333	+	new_left,
1334	+	new_cigars,
1335	+	bh.antisense_align(),
1336	+	bh.antisense_splice(),
1337	+	0, /* edit_dist - needs to be recalculated */
1338	+	0, /* splice_mms - needs to be recalculated */
1339	+	false);
1340	+
1341	+	new_bh.seq(bh.seq());
1342	+	new_bh.qual(bh.qual());
1343	+
1344	+	vector<string> aux_fields;
1345	+	bowtie_sam_extra(new_bh, rt, aux_fields);
1346	+
1347	+	vector<string>::const_iterator aux_iter = aux_fields.begin();
1348	+	for (; aux_iter != aux_fields.end(); ++aux_iter)
1349	+	{
1350	+	const string& aux_field = *aux_iter;
1351	+	if (strncmp(aux_field.c_str(), "AS", 2) == 0)
1352	+	{
1353	+	int alignment_score = atoi(aux_field.c_str() + 5);
1354	+	new_bh.alignment_score(alignment_score);
1355	+	}
1356	+	else if (strncmp(aux_field.c_str(), "XM", 2) == 0)
1357	+	{
1358	+	int XM_value = atoi(aux_field.c_str() + 5);
1359	+	new_bh.edit_dist(XM_value);
1360	+	}
1361	+	}
1362	+
1363	+	/*
1364	+	fprintf(stderr, "\t%d %d-%d %s (AS:%d XM:%d) with junc %u-%u\n",
1365	+	new_bh.insert_id(), new_bh.left(), new_bh.right(),
1366	+	print_cigar(new_bh.cigar()).c_str(),
1367	+	new_bh.alignment_score(), new_bh.edit_dist(),
1368	+	junc.left, junc.right);
1369	+	*/
1370	+
1371	+	++lb;
1372	+	}
1373	+	}
1374	+	#endif
1375	+	}
1376	+	}
1377	+
1378	+	switch(op.opcode)
1379	+	{
1380	+	case REF_SKIP:
1381	+	pos += op.length;
1382	+	break;
1383	+	case rEF_SKIP:
1384	+	pos -= op.length;
1385	+	break;
1386	+	case MATCH:
1387	+	case DEL:
1388	+	pos += op.length;
1389	+	break;
1390	+	case mATCH:
1391	+	case dEL:
1392	+	pos -= op.length;
1393	+	break;
1394	+	case FUSION_FF:
1395	+	case FUSION_FR:
1396	+	case FUSION_RF:
1397	+	case FUSION_RR:
1398	+	pos = op.length;
1399	+	refid = bh.ref_id2();
1400	+	break;
1401	+	default:
1402	+	break;
1403	+	}
1404	+	}
1405	+	}
1406	+
1407	+	hits.hits.insert(hits.hits.end(), additional_hits.begin(), additional_hits.end());
1408	+
1409	+	std::sort(hits.hits.begin(), hits.hits.end());
1410	+	vector<BowtieHit>::iterator new_end = std::unique(hits.hits.begin(), hits.hits.end());
1411	+	hits.hits.erase(new_end, hits.hits.end());
1412	+	}
1413	+
1414	+
1415		// events include splice junction, indels, and fusion points.
1416		struct ConsensusEventsWorker
1417		{
#	Line 1042 \| Line 1556
1556		{
1557		void operator()()
1558		{
1559	+	rng.seed(1);
1560	+
1561		ReadTable it;
1562		GBamWriter bam_writer(bam_output_fname.c_str(), sam_header.c_str());
1563
#	Line 1059 \| Line 1575
1575		GBamWriter* left_um_out=new GBamWriter(left_um_fname.c_str(), sam_header.c_str());
1576		GBamWriter* right_um_out=NULL;
1577
1062	–	//if (left_um_out.openWrite(left_um_fname.c_str(), zpacker)==NULL)
1063	–	// err_die("Error: cannot open file %s for writing!\n",left_um_fname.c_str());
1578		ReadStream right_reads_file(right_reads_fname);
1579		if (right_reads_offset > 0)
1580		right_reads_file.seek(right_reads_offset);
#	Line 1113 \| Line 1627
1627		{
1628		HitsForRead best_hits;
1629		best_hits.insert_id = curr_left_obs_order;
1630	<	bool got_read=left_reads_file.getRead(curr_left_obs_order, l_read,
1631	<	reads_format, false, begin_id, end_id,
1632	<	left_um_out, false);
1633	<	//assert(got_read);
1120	<
1121	<	if (right_reads_file.file()) {
1122	<	got_read=right_reads_file.getRead(curr_left_obs_order, r_read,
1123	<	reads_format, false, begin_id, end_id,
1124	<	right_um_out, true);
1125	<	//assert(got_read);
1126	<	}
1127	<
1630	>
1631	>	realign_reads(curr_left_hit_group, rt, junctions, *rev_junctions,
1632	>	insertions, deletions, rev_deletions, fusions);
1633	>
1634		exclude_hits_on_filtered_junctions(*junctions, curr_left_hit_group);
1635
1636		// Process hits for left singleton, select best alignments
#	Line 1134 \| Line 1640
1640
1641		if (best_hits.hits.size() > 0)
1642		{
1643	<	update_junctions(best_hits, *final_junctions);
1644	<	update_insertions_and_deletions(best_hits, final_insertions, final_deletions);
1645	<	update_fusions(best_hits, rt, final_fusions, *fusions);
1646	<
1647	<	print_sam_for_single(*rt,
1648	<	best_hits,
1649	<	(right_map_fname.empty() ? FRAG_UNPAIRED : FRAG_LEFT),
1650	<	l_read,
1651	<	bam_writer,
1652	<	left_um_out);
1643	>	bool got_read = left_reads_file.getRead(curr_left_obs_order, l_read,
1644	>	reads_format, false, begin_id, end_id,
1645	>	left_um_out, false);
1646	>	assert (got_read);
1647	>
1648	>	if (got_read)
1649	>	{
1650	>	update_junctions(best_hits, *final_junctions);
1651	>	update_insertions_and_deletions(best_hits, final_insertions, final_deletions);
1652	>	update_fusions(best_hits, rt, final_fusions, *fusions);
1653	>
1654	>	print_sam_for_single(*rt,
1655	>	best_hits,
1656	>	(right_map_fname.empty() ? FRAG_UNPAIRED : FRAG_LEFT),
1657	>	l_read.alt_name,
1658	>	bam_writer,
1659	>	rng);
1660	>	}
1661		}
1662	+
1663		// Get next hit group
1664		left_hs.next_read_hits(curr_left_hit_group);
1665		curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id);
#	Line 1160 \| Line 1675
1675
1676		if (curr_right_obs_order >= begin_id)
1677		{
1678	<	bool got_read=right_reads_file.getRead(curr_right_obs_order, r_read,
1679	<	reads_format, false, begin_id, end_id,
1165	<	right_um_out, false);
1166	<
1167	<	//assert(got_read);
1168	<	/*
1169	<	fprintf(right_um_out.file, "@%s #MAPPED#\n%s\n+\n%s\n", r_read.alt_name.c_str(),
1170	<	r_read.seq.c_str(), r_read.qual.c_str());
1171	<	*/
1172	<	got_read=left_reads_file.getRead(curr_right_obs_order, l_read,
1173	<	reads_format, false, begin_id, end_id,
1174	<	left_um_out, true);
1175	<	//assert(got_read);
1176	<
1678	>	realign_reads(curr_right_hit_group, rt, junctions, *rev_junctions,
1679	>	insertions, deletions, rev_deletions, fusions);
1680		exclude_hits_on_filtered_junctions(*junctions, curr_right_hit_group);
1681
1682		// Process hit for right singleton, select best alignments
#	Line 1183 \| Line 1686
1686
1687		if (best_hits.hits.size() > 0)
1688		{
1689	<	update_junctions(best_hits, *final_junctions);
1690	<	update_insertions_and_deletions(best_hits, final_insertions, final_deletions);
1691	<	update_fusions(best_hits, rt, final_fusions, *fusions);
1692	<
1693	<	print_sam_for_single(*rt,
1694	<	best_hits,
1695	<	FRAG_RIGHT,
1696	<	r_read,
1697	<	bam_writer,
1698	<	right_um_out);
1689	>	bool got_read =right_reads_file.getRead(curr_right_obs_order, r_read,
1690	>	reads_format, false, begin_id, end_id,
1691	>	right_um_out, false);
1692	>	assert (got_read);
1693	>
1694	>	if (got_read)
1695	>	{
1696	>	update_junctions(best_hits, *final_junctions);
1697	>	update_insertions_and_deletions(best_hits, final_insertions, final_deletions);
1698	>	update_fusions(best_hits, rt, final_fusions, *fusions);
1699	>
1700	>	print_sam_for_single(*rt,
1701	>	best_hits,
1702	>	FRAG_RIGHT,
1703	>	r_read.alt_name,
1704	>	bam_writer,
1705	>	rng);
1706	>	}
1707		}
1708		}
1709
#	Line 1207 \| Line 1718
1718		curr_left_obs_order < end_id &&
1719		curr_left_obs_order != VMAXINT32)
1720		{
1721	+	realign_reads(curr_left_hit_group, rt, junctions, *rev_junctions,
1722	+	insertions, deletions, rev_deletions, fusions);
1723		exclude_hits_on_filtered_junctions(*junctions, curr_left_hit_group);
1724	+
1725	+	realign_reads(curr_right_hit_group, rt, junctions, *rev_junctions,
1726	+	insertions, deletions, rev_deletions, fusions);
1727		exclude_hits_on_filtered_junctions(*junctions, curr_right_hit_group);
1728	+
1729		if (curr_left_hit_group.hits.empty())
1730		{ //only right read mapped
1731		//write it in the mapped file with the #MAPPED# flag
1215	–
1216	–	bool got_read=right_reads_file.getRead(curr_left_obs_order, r_read,
1217	–	reads_format, false, begin_id, end_id,
1218	–	right_um_out, false);
1219	–	//assert(got_read);
1220	–	/*
1221	–	fprintf(right_um_out.file, "@%s #MAPPED#\n%s\n+\n%s\n", r_read.alt_name.c_str(),
1222	–	r_read.seq.c_str(), r_read.qual.c_str());
1223	–	*/
1732		HitsForRead right_best_hits;
1733		right_best_hits.insert_id = curr_right_obs_order;
1734
#	Line 1230 \| Line 1738
1738
1739		if (right_best_hits.hits.size() > 0)
1740		{
1741	<	update_junctions(right_best_hits, *final_junctions);
1742	<	update_insertions_and_deletions(right_best_hits, final_insertions, final_deletions);
1743	<	update_fusions(right_best_hits, rt, final_fusions, *fusions);
1744	<
1745	<	print_sam_for_single(*rt,
1746	<	right_best_hits,
1747	<	FRAG_RIGHT,
1748	<	r_read,
1749	<	bam_writer,
1750	<	right_um_out);
1741	>	bool got_read = right_reads_file.getRead(curr_left_obs_order, r_read,
1742	>	reads_format, false, begin_id, end_id,
1743	>	right_um_out, false);
1744	>	assert (got_read);
1745	>
1746	>	if (got_read)
1747	>	{
1748	>	update_junctions(right_best_hits, *final_junctions);
1749	>	update_insertions_and_deletions(right_best_hits, final_insertions, final_deletions);
1750	>	update_fusions(right_best_hits, rt, final_fusions, *fusions);
1751	>
1752	>	print_sam_for_single(*rt,
1753	>	right_best_hits,
1754	>	FRAG_RIGHT,
1755	>	r_read.alt_name,
1756	>	bam_writer,
1757	>	rng);
1758	>	}
1759		}
1760		}
1761		else if (curr_right_hit_group.hits.empty())
#	Line 1247 \| Line 1763
1763		HitsForRead left_best_hits;
1764		left_best_hits.insert_id = curr_left_obs_order;
1765		//only left read mapped
1250	–	bool got_read=left_reads_file.getRead(curr_left_obs_order, l_read,
1251	–	reads_format, false, begin_id, end_id,
1252	–	left_um_out, false);
1253	–
1766		// Process hits for left singleton, select best alignments
1767		read_best_alignments(curr_left_hit_group, left_best_hits, *gtf_junctions,
1768		junctions, insertions, deletions, fusions, *coverage,
#	Line 1258 \| Line 1770
1770
1771		if (left_best_hits.hits.size() > 0)
1772		{
1773	<	update_junctions(left_best_hits, *final_junctions);
1774	<	update_insertions_and_deletions(left_best_hits, final_insertions, final_deletions);
1775	<	update_fusions(left_best_hits, rt, final_fusions, *fusions);
1776	<
1777	<	print_sam_for_single(*rt,
1778	<	left_best_hits,
1779	<	FRAG_LEFT,
1780	<	l_read,
1781	<	bam_writer,
1782	<	left_um_out);
1773	>	bool got_read =left_reads_file.getRead(curr_left_obs_order, l_read,
1774	>	reads_format, false, begin_id, end_id,
1775	>	left_um_out, false);
1776	>	assert (got_read);
1777	>
1778	>	if (got_read)
1779	>	{
1780	>	update_junctions(left_best_hits, *final_junctions);
1781	>	update_insertions_and_deletions(left_best_hits, final_insertions, final_deletions);
1782	>	update_fusions(left_best_hits, rt, final_fusions, *fusions);
1783	>
1784	>	print_sam_for_single(*rt,
1785	>	left_best_hits,
1786	>	FRAG_LEFT,
1787	>	l_read.alt_name,
1788	>	bam_writer,
1789	>	rng);
1790	>	}
1791		}
1792		}
1793		else
#	Line 1292 \| Line 1812
1812
1813		if (best_hits.size() > 0)
1814		{
1815	<	update_junctions(best_left_hit_group, *final_junctions);
1816	<	update_insertions_and_deletions(best_left_hit_group, final_insertions, final_deletions);
1817	<	update_fusions(best_left_hit_group, rt, final_fusions, *fusions);
1818	<
1819	<	update_junctions(best_right_hit_group, *final_junctions);
1820	<	update_insertions_and_deletions(best_right_hit_group, final_insertions, final_deletions);
1821	<	update_fusions(best_right_hit_group, rt, final_fusions, *fusions);
1822	<
1823	<	pair_support(best_hits, final_fusions, fusions);
1824	<
1825	<	print_sam_for_pair(*rt,
1826	<	best_hits,
1827	<	grade,
1828	<	left_reads_file,
1829	<	right_reads_file,
1830	<	bam_writer,
1831	<	left_um_out,
1832	<	right_um_out,
1833	<	begin_id,
1834	<	end_id);
1815	>	bool got_left_read = left_reads_file.getRead(best_hits[0].first.insert_id(), l_read,
1816	>	reads_format, false, begin_id, end_id,
1817	>	left_um_out, false);
1818	>
1819	>	bool got_right_read = right_reads_file.getRead(best_hits[0].first.insert_id(), r_read,
1820	>	reads_format, false, begin_id, end_id,
1821	>	right_um_out, false);
1822	>
1823	>	assert (got_left_read && got_right_read);
1824	>
1825	>	if (got_left_read && got_right_read)
1826	>	{
1827	>	update_junctions(best_left_hit_group, *final_junctions);
1828	>	update_insertions_and_deletions(best_left_hit_group, final_insertions, final_deletions);
1829	>	update_fusions(best_left_hit_group, rt, final_fusions, *fusions);
1830	>
1831	>	update_junctions(best_right_hit_group, *final_junctions);
1832	>	update_insertions_and_deletions(best_right_hit_group, final_insertions, final_deletions);
1833	>	update_fusions(best_right_hit_group, rt, final_fusions, *fusions);
1834	>
1835	>	pair_support(best_hits, final_fusions, fusions);
1836	>
1837	>	print_sam_for_pair(*rt,
1838	>	best_hits,
1839	>	grade,
1840	>	bam_writer,
1841	>	l_read.alt_name,
1842	>	r_read.alt_name,
1843	>	rng,
1844	>	begin_id,
1845	>	end_id);
1846	>	}
1847		}
1848		}
1849
#	Line 1359 \| Line 1891
1891		JunctionSet* gtf_junctions;
1892
1893		JunctionSet* junctions;
1894	+	JunctionSet* rev_junctions;
1895		InsertionSet* insertions;
1896		DeletionSet* deletions;
1897	+	DeletionSet* rev_deletions;
1898		FusionSet* fusions;
1899		Coverage* coverage;
1900
#	Line 1377 \| Line 1911
1911		int64_t left_map_offset;
1912		int64_t right_reads_offset;
1913		int64_t right_map_offset;
1914	+
1915	+	boost::random::mt19937 rng;
1916		};
1917
1918		void driver(const string& bam_output_fname,
#	Line 1394 \| Line 1930
1930		num_threads = 1;
1931
1932		RefSequenceTable rt(sam_header, true);
1933	<
1398	<	if (fusion_search)
1399	<	get_seqs(ref_stream, rt, true);
1933	>	get_seqs(ref_stream, rt, true);
1934
1935		srandom(1);
1936
#	Line 1513 \| Line 2047
2047		DeletionSet& deletions = vdeletions[0];
2048		FusionSet& fusions = vfusions[0];
2049		Coverage& coverage = vcoverages[0];
2050	<	for (size_t i = 1; i < num_threads; ++i)
2050	>	for (int i = 1; i < num_threads; ++i)
2051		{
2052		merge_with(junctions, vjunctions[i]);
2053		vjunctions[i].clear();
#	Line 1532 \| Line 2066
2066		}
2067
2068		coverage.calculate_coverage();
2069	+
2070	+	JunctionSet rev_junctions;
2071	+	JunctionSet::const_iterator junction_iter = junctions.begin();
2072	+	for (; junction_iter != junctions.end(); ++junction_iter)
2073	+	{
2074	+	const Junction& junction = junction_iter->first;
2075	+	Junction rev_junction = junction;
2076	+	rev_junction.left = junction.right;
2077	+	rev_junction.right = junction.left;
2078	+	rev_junctions[rev_junction] = junction_iter->second;
2079	+	}
2080	+
2081	+	DeletionSet rev_deletions;
2082	+	#if 0
2083	+	DeletionSet::const_iterator deletion_iter = deletions.begin();
2084	+	for (; deletion_iter != deletions.end(); ++deletion_iter)
2085	+	{
2086	+	const Deletion& deletion = deletion_iter->first;
2087	+	Deletion rev_deletion = deletion;
2088	+	rev_deletion.left = deletion.right;
2089	+	rev_deletion.right = deletion.left;
2090	+	rev_deletions[rev_deletion] = deletion_iter->second;
2091	+	}
2092	+	#endif
2093
2094		size_t num_unfiltered_juncs = junctions.size();
2095		fprintf(stderr, "Loaded %lu junctions\n", (long unsigned int) num_unfiltered_juncs);
#	Line 1590 \| Line 2148
2148
2149		worker.gtf_junctions = &gtf_junctions;
2150		worker.junctions = &junctions;
2151	+	worker.rev_junctions = &rev_junctions;
2152		worker.insertions = &insertions;
2153		worker.deletions = &deletions;
2154	+	worker.rev_deletions = &rev_deletions;
2155		worker.fusions = &fusions;
2156		worker.coverage = &coverage;
2157
#	Line 1644 \| Line 2204
2204		InsertionSet final_insertions = vfinal_insertions[0];
2205		DeletionSet final_deletions = vfinal_deletions[0];
2206		FusionSet final_fusions = vfinal_fusions[0];
2207	<	for (size_t i = 1; i < num_threads; ++i)
2207	>	for (int i = 1; i < num_threads; ++i)
2208		{
2209		merge_with(final_junctions, vfinal_junctions[i]);
2210		vfinal_junctions[i].clear();
#	Line 1780 \| Line 2340
2340
2341		string left_reads_filename = argv[optind++];
2342		string unzcmd=getUnpackCmd(left_reads_filename, false);
1783	–
2343		string right_map_filename;
2344		string right_reads_filename;
2345

Diff Legend

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