ViewVC Help

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

root/gclib/tophat_cpp/bwt_map.cpp

#	Line 373 | Line 373
373
374		int anchor_mismatch = 0;
375
376	+
377	+	void parseSegReadName(char* name, char*& name_tags, bool strip_slash,
378	+	bool &end, unsigned int &seg_offset, unsigned int& seg_num,
379	+	unsigned int & num_segs) {
380	+	char* pipe = strrchr(name, '|');
381	+	if (pipe)
382	+	{
383	+	if (name_tags)
384	+	strcpy(name_tags, pipe);
385	+
386	+	char* tag_buf = pipe + 1;
387	+	if (strchr(tag_buf, ':'))
388	+	{
389	+	sscanf(tag_buf, "%u:%u:%u", &seg_offset, &seg_num, &num_segs);
390	+	if (seg_num + 1 == num_segs)
391	+	end = true;
392	+	else
393	+	end = false;
394	+	}
395	+
396	+	*pipe = 0;
397	+	}
398	+
399	+	// Stripping the slash and number following it gives the insert name
400	+	char* slash = strrchr(name, '/');
401	+	if (strip_slash && slash)
402	+	*slash = 0;
403	+	}
404	+
405	+
406		bool SplicedBowtieHitFactory::get_hit_from_buf(const char* orig_bwt_buf,
407		BowtieHit& bh,
408		bool strip_slash,
#	Line 426 | Line 456
456		unsigned int num_segs = 0;
457
458		// Copy the tag out of the name field before we might wipe it out
459	<	char* pipe = strrchr(name, '|');
430	<	if (pipe)
431	<	{
432	<	if (name_tags)
433	<	strcpy(name_tags, pipe);
459	>	parseSegReadName(name, name_tags, strip_slash, end, seg_offset, seg_num, num_segs);
460
435	–	char* tag_buf = pipe + 1;
436	–	if (strchr(tag_buf, ':'))
437	–	{
438	–	sscanf(tag_buf, "%u:%u:%u", &seg_offset, &seg_num, &num_segs);
439	–	if (seg_num + 1 == num_segs)
440	–	end = true;
441	–	else
442	–	end = false;
443	–	}
444	–
445	–	*pipe = 0;
446	–	}
447	–	// Stripping the slash and number following it gives the insert name
448	–	char* slash = strrchr(name, '/');
449	–	if (strip_slash && slash)
450	–	*slash = 0;
451	–
452	–	//int read_len = strlen(sequence);
453	–
461		// Add this alignment to the table of hits for this half of the
462		// Bowtie map
463
#	Line 712 | Line 719
719		return false;
720		}
721
722	+
723	+	char parseCigar(vector<CigarOp>& cigar, const char* cigar_str,
724	+	bool &spliced_alignment)   {
725	+	// Mostly pilfered direct from the SAM tools:
726	+	const char* p_cig = cigar_str;
727	+	while (*p_cig)
728	+	{
729	+	char* t;
730	+	int length = (int)strtol(p_cig, &t, 10);
731	+	if (length <= 0)
732	+	{
733	+	//fprintf (stderr, "CIGAR op has zero length\n");
734	+	return false;
735	+	}
736	+	char op_char = toupper(*t);
737	+	CigarOpCode opcode;
738	+	if (op_char == 'M') opcode = MATCH;
739	+	else if (op_char == 'I') opcode = INS;
740	+	else if (op_char == 'D') opcode = DEL;
741	+	else if (op_char == 'N')
742	+	{
743	+	if (length > max_report_intron_length)
744	+	return false;
745	+	opcode = REF_SKIP;
746	+	spliced_alignment = true;
747	+	}
748	+	else if (op_char == 'S') opcode = SOFT_CLIP;
749	+	else if (op_char == 'H') opcode = HARD_CLIP;
750	+	else if (op_char == 'P') opcode = PAD;
751	+	else
752	+	{
753	+	fprintf (stderr, "invalid CIGAR operation\n");
754	+	return false;
755	+	}
756	+	p_cig = t + 1;
757	+	//i += length;
758	+	cigar.push_back(CigarOp(opcode, length));
759	+	}
760	+	if (*p_cig) {
761	+	fprintf (stderr, "unmatched CIGAR operation\n");
762	+	return *p_cig;
763	+	}
764	+	return 0;
765	+	}
766	+
767	+	void getSAMmismatches(char* &buf, vector<CigarOp>& cigar,
768	+	int* mismatches, int& num_mismatches, int& sam_nm, bool& antisense_splice) {
769	+	const char* tag_buf = buf;
770	+	sam_nm=0;
771	+	num_mismatches=0;
772	+	while((tag_buf = get_token((char**)&buf,"\t")))
773	+	{
774	+	vector<string> tuple_fields;
775	+	tokenize(tag_buf,":", tuple_fields);
776	+	if (tuple_fields.size() == 3)
777	+	{
778	+	if (tuple_fields[0] == "XS")
779	+	{
780	+	if (tuple_fields[2] == "-")
781	+	antisense_splice = true;
782	+	}
783	+	else if (tuple_fields[0] == "NM")
784	+	{
785	+	sam_nm = atoi(tuple_fields[2].c_str());
786	+	}
787	+	else if (tuple_fields[0] == "NS")
788	+	{
789	+	//ignored for now
790	+	}
791	+	else if (tuple_fields[0] == "MD")
792	+	{
793	+	const char* p=tuple_fields[2].c_str();
794	+	int bi=0; //base offset position in the read
795	+	while (*p != 0) {
796	+	if (isdigit(*p)) {
797	+	int v=atoi(p);
798	+	do { p++; } while (isdigit(*p));
799	+	bi+=v;
800	+	}
801	+	while (isupper(*p)) {
802	+	p++;
803	+	mismatches[num_mismatches]=bi;
804	+	num_mismatches++;
805	+	bi++;
806	+	}
807	+	if (*p=='^') { //reference deletion
808	+	p++;
809	+	while (isupper(*p)) { //insert read bases
810	+	p++; bi++;
811	+	}
812	+	}
813	+	}
814	+	}
815	+	else
816	+	{
817	+	//fprintf(stderr, "%s attribute not supported\n", tuple_fields[0].c_str());
818	+	//return false;
819	+	}
820	+	}
821	+	}
822	+	/* By convention,the NM field of the SAM record
823	+	*  counts an insertion or deletion. I dont' think
824	+	*  we want the mismatch count in the BowtieHit
825	+	*  record to reflect this. Therefore, subtract out
826	+	*  the mismatches due to in/dels
827	+	*/
828	+	for(vector<CigarOp>::const_iterator itr = cigar.begin(); itr != cigar.end(); ++itr){
829	+	if(itr->opcode == INS){
830	+	sam_nm -= itr->length;
831	+	}
832	+	if(itr->opcode == DEL){
833	+	sam_nm -= itr->length;
834	+	}
835	+	}
836	+	}
837	+
838		bool SAMHitFactory::get_hit_from_buf(const char* orig_bwt_buf,
839		BowtieHit& bh,
840		bool strip_slash,
#	Line 773 | Line 896
896		unsigned int num_segs = 0;
897
898		// Copy the tag out of the name field before we might wipe it out
899	<	char* pipe = strrchr(name, '|');
777	<	if (pipe)
778	<	{
779	<	if (name_tags)
780	<	strcpy(name_tags, pipe);
781	<
782	<	char* tag_buf = pipe + 1;
783	<	if (strchr(tag_buf, ':'))
784	<	{
785	<	sscanf(tag_buf, "%u:%u:%u", &seg_offset, &seg_num, &num_segs);
786	<	if (seg_num + 1 == num_segs)
787	<	end = true;
788	<	else
789	<	end = false;
790	<	}
791	<
792	<	*pipe = 0;
793	<	}
794	<
795	<	// Stripping the slash and number following it gives the insert name
796	<	char* slash = strrchr(name, '/');
797	<	if (strip_slash && slash)
798	<	*slash = 0;
899	>	parseSegReadName(name, name_tags, strip_slash, end, seg_offset, seg_num, num_segs);
900
800	–	char* p_cig = cigar_str;
801	–	//int len = strlen(sequence);
901		vector<CigarOp> cigar;
902		bool spliced_alignment = false;
903	<	// Mostly pilfered direct from the SAM tools:
904	<	while (*p_cig)
905	<	{
807	<	char* t;
808	<	int length = (int)strtol(p_cig, &t, 10);
809	<	if (length <= 0)
810	<	{
811	<	//fprintf (stderr, "CIGAR op has zero length\n");
812	<	return false;
813	<	}
814	<	char op_char = toupper(*t);
815	<	CigarOpCode opcode;
816	<	if (op_char == 'M') opcode = MATCH;
817	<	else if (op_char == 'I') opcode = INS;
818	<	else if (op_char == 'D') opcode = DEL;
819	<	else if (op_char == 'N')
820	<	{
821	<	if (length > max_report_intron_length)
822	<	return false;
823	<	opcode = REF_SKIP;
824	<	spliced_alignment = true;
825	<	}
826	<	else if (op_char == 'S') opcode = SOFT_CLIP;
827	<	else if (op_char == 'H') opcode = HARD_CLIP;
828	<	else if (op_char == 'P') opcode = PAD;
829	<	else
830	<	{
831	<	fprintf (stderr, "invalid CIGAR operation\n");
832	<	return false;
833	<	}
834	<	p_cig = t + 1;
835	<	//i += length;
836	<	cigar.push_back(CigarOp(opcode, length));
837	<	}
838	<	if (*p_cig)
839	<	{
840	<	fprintf (stderr, "unmatched CIGAR operation\n");
841	<	return false;
842	<	}
903	>	if (parseCigar(cigar, cigar_str, spliced_alignment)) {
904	>	return false;
905	>	}
906
907		//vector<string> attributes;
908		//tokenize(tag_buf, " \t",attributes);
909	<
909	>
910		bool antisense_splice = false;
911	<	unsigned char num_mismatches = 0;
912	<	unsigned char num_splice_anchor_mismatches = 0;
913	<	const char* tag_buf = buf;
914	<
852	<	//      while((tag_buf = get_token((char**)&buf,"\t")))
853	<	//      {
854	<	//              vector<string> tuple_fields;
855	<	//              tokenize(tag_buf,":", tuple_fields);
856	<	//              if (tuple_fields.size() == 3)
857	<	//              {
858	<	//                      if (tuple_fields[0] == "XS")
859	<	//                      {
860	<	//                              if (tuple_fields[2] == "-")
861	<	//                                      antisense_splice = true;
862	<	//                      }
863	<	//                      else if (tuple_fields[0] == "NM")
864	<	//                      {
865	<	//                              num_mismatches = atoi(tuple_fields[2].c_str());
866	<	//                      }
867	<	//                      else if (tuple_fields[0] == "NS")
868	<	//                      {
869	<	//                              num_splice_anchor_mismatches = atoi(tuple_fields[2].c_str());
870	<	//                      }
871	<	//                      else
872	<	//                      {
873	<	//                              fprintf(stderr, "%s attribute not supported\n", tuple_fields[0].c_str());
874	<	//                              return false;
875	<	//                      }
876	<	//              }
877	<	//      }
878	<
879	<	while((tag_buf = get_token((char**)&buf,"\t")))
880	<	{
881	<	vector<string> tuple_fields;
882	<	tokenize(tag_buf,":", tuple_fields);
883	<	if (tuple_fields.size() == 3)
884	<	{
885	<	if (tuple_fields[0] == "XS")
886	<	{
887	<	if (tuple_fields[2] == "-")
888	<	antisense_splice = true;
889	<	}
890	<	else if (tuple_fields[0] == "NM")
891	<	{
892	<	num_mismatches = atoi(tuple_fields[2].c_str());
893	<	}
894	<	else if (tuple_fields[0] == "NS")
895	<	{
896	<	//ignored for now
897	<	}
898	<	else
899	<	{
900	<	//fprintf(stderr, "%s attribute not supported\n", tuple_fields[0].c_str());
901	<	//return false;
902	<	}
903	<	}
904	<	}
911	>	int num_mismatches = 0;
912	>	int sam_nm = 0; //the value of the NM tag (edit distance)
913	>	int mismatches[1024];
914	>	getSAMmismatches(buf, cigar, mismatches, num_mismatches, sam_nm, antisense_splice);
915
906	–
907	–	/*
908	–	* By convention,the NM field of the SAM record
909	–	* counts an insertion or deletion. I dont' think
910	–	* we want the mismatch count in the BowtieHit
911	–	* record to reflect this. Therefore, subtract out
912	–	* the mismatches due to in/dels
913	–	*/
914	–	for(vector<CigarOp>::const_iterator itr = cigar.begin(); itr != cigar.end(); ++itr){
915	–	if(itr->opcode == INS){
916	–	num_mismatches -= itr->length;
917	–	}
918	–	if(itr->opcode == DEL){
919	–	num_mismatches -= itr->length;
920	–	}
921	–	}
922	–	//      vector<string> toks;
923	–	//      tokenize(tag_buf, ":", toks);
916		if (spliced_alignment)
917		{
918		bh = create_hit(name,
#	Line 930 | Line 922
922		sam_flag & 0x0010,
923		antisense_splice,
924		num_mismatches,
925	<	num_splice_anchor_mismatches,
925	>	0,
926		end);
935	–	return true;
936	–
927		}
928		else
929		{
930		//assert(cigar.size() == 1 && cigar[0].opcode == MATCH);
941	–
931		bh = create_hit(name,
932		text_name,
933		text_offset - 1, // SAM files are 1-indexed
#	Line 948 | Line 937
937		num_mismatches,
938		0,
939		end);
951	–	return true;
940		}
941	<
942	<	return false;
941	>	return true;
942	>	}
943	>
944	>
945	>	void spliceCigar(vector<CigarOp>& splcigar, vector<CigarOp>& cigar,
946	>	int8_t l_len, int mlen, int8_t r_len, CigarOpCode opcode) {
947	>	//merge the original 'cigar' with the new insert/gap operation
948	>	// at position l_len into splcigar;
949	>	//find the original opcode that croses l_len location
950	>	int ofs=0;
951	>	bool found=false;
952	>	for (size_t c = 0 ; c < cigar.size(); ++c)
953	>	{
954	>	ofs+=cigar[c].length;
955	>	if (!found) {
956	>	if (ofs>=l_len) {
957	>	found=true;
958	>	//inject new code here
959	>	}
960	>	else {
961	>	//not found yet, just copy these codes
962	>	splcigar.push_back(cigar[c]);
963	>	}
964	>	}
965	>	else { //found yet
966	>
967	>	}
968	>	}
969		}
970
971		bool SplicedSAMHitFactory::get_hit_from_buf(const char* orig_bwt_buf,
#	Line 971 | Line 985
985		if (bwt_buf[0] == '@')
986		return false;
987
974	–	//char orientation;
975	–	/*
976	–	char mismatches[buf_size];
977	–
978	–	//char* orientation_str = get_token((char**)&buf,"\t");
979	–
980	–	mismatches[0] = 0;
981	–	char* mismatches_str = get_token((char**)&buf, "\t");
982	–	if (mismatches_str)
983	–	strcpy(mismatches, mismatches_str);
984	–
985	–	orientation = orientation_str[0];
986	–	text_offset = atoi(text_offset_str);
987	–	*/
988	–
988		char* buf = bwt_buf;
989		char* name = get_token((char**)&buf,"\t");
990		char* sam_flag_str = get_token((char**)&buf,"\t");
#	Line 996 | Line 995
995		const char* mate_ref_str =  get_token((char**)&buf,"\t");
996		const char* mate_pos_str =  get_token((char**)&buf,"\t");
997		const char* inferred_insert_sz_str =  get_token((char**)&buf,"\t");
998	<
998	>	int num_mismatches=0;
999	>	int mismatches[1024]; //list of 0-based mismatch positions in this read
1000	>	//parsed from SAM's MD:Z: tag
1001		const char* seq_str =  get_token((char**)&buf,"\t");
1002		if (seq)
1003		strcpy(seq, seq_str);
#	Line 1024 | Line 1025
1025
1026		int sam_flag = atoi(sam_flag_str);
1027		int text_offset = atoi(text_offset_str);
1028	<
1028	<	//##############################################
1028	>	text_offset--; //SAM is 1-based, Bowtie is 0-based
1029		bool end = true;
1030		unsigned int seg_offset = 0;
1031		unsigned int seg_num = 0;
1032		unsigned int num_segs = 0;
1033
1034		// Copy the tag out of the name field before we might wipe it out
1035	<	char* pipe = strrchr(name, '|');
1036	<	if (pipe)
1037	<	{
1038	<	if (name_tags)
1039	<	strcpy(name_tags, pipe);
1035	>	parseSegReadName(name, name_tags, strip_slash, end, seg_offset, seg_num, num_segs);
1036
1037	<	char* tag_buf = pipe + 1;
1038	<	if (strchr(tag_buf, ':'))
1039	<	{
1040	<	sscanf(tag_buf, "%u:%u:%u", &seg_offset, &seg_num, &num_segs);
1041	<	if (seg_num + 1 == num_segs)
1046	<	end = true;
1047	<	else
1048	<	end = false;
1049	<	}
1037	>	vector<CigarOp> samcigar;
1038	>	bool spliced_alignment = false;
1039	>	if (parseCigar(samcigar, cigar_str, spliced_alignment)) {
1040	>	return false;
1041	>	}
1042
1043	<	*pipe = 0;
1044	<	}
1045	<	// Stripping the slash and number following it gives the insert name
1054	<	char* slash = strrchr(name, '/');
1055	<	if (strip_slash && slash)
1056	<	*slash = 0;
1043	>	bool antisense_splice = false;
1044	>	int sam_nm = 0;
1045	>	getSAMmismatches(buf, samcigar, mismatches, num_mismatches, sam_nm, antisense_splice);
1046
1047	<	//int read_len = strlen(sequence);
1047	>	//##############################################
1048
1049		// Add this alignment to the table of hits for this half of the
1050		// Bowtie map
#	Line 1087 | Line 1076
1076		if (splice_toks.size() != 2)
1077		{
1078		fprintf(stderr, "Warning: found malformed splice record, skipping:\n");
1079	<	//fprintf(stderr, "%s (token: %s)\n", text_name,
1079	>	//fprintf(stderr, "\t%s (token: %s)\n", text_name,
1080		//        toks[num_extra_toks + splice_field].c_str());
1081		return false;
1082		}
1083	<
1083	>	int spl_num_mismatches=0;
1084		//
1085		// check for an insertion hit
1086		//
1087		if(toks[num_extra_toks + junction_type_field] == "ins")
1088		{
1089		int8_t spliced_read_len = strlen(seq_str);
1090	<	/*
1091	<	* The 0-based position of the left edge of the alignment. Note that this
1092	<	* value may need to be futher corrected to account for the presence of
1104	<	* of the insertion.
1105	<	*/
1090	>	// The 0-based position of the left edge of the alignment. Note that this
1091	>	// value may need to be futher corrected to account for the presence of
1092	>	// of the insertion.
1093		uint32_t left = atoi(toks[num_extra_toks + left_window_edge_field].c_str()) + text_offset;
1094		uint32_t right = left + spliced_read_len - 1;
1095
1096	<	/*
1110	<	* The 0-based position of the last genomic sequence before the insertion
1111	<	*/
1096	>	// The 0-based position of the last genomic sequence before the insertion
1097		uint32_t left_splice_pos = atoi(splice_toks[0].c_str());
1098
1099		string insertedSequence = splice_toks[1];
1100	<	/*
1101	<	* The 0-based position of the first genomic sequence after teh insertion
1117	<	*/
1100	>	// The 0-based position of the first genomic sequence after the insertion
1101	>
1102		uint32_t right_splice_pos = left_splice_pos + 1;
1103		if(left > left_splice_pos){
1104		/*
#	Line 1157 | Line 1141
1141		return false;
1142		}
1143
1144	<	char* pch = strtok( mismatches, ",");
1145	<	unsigned char num_mismatches = 0;
1144	>	//char* pch = strtok( mismatches, ",");
1145	>	//unsigned char num_mismatches = 0;
1146	>	//text_offset holds the left end of the alignment,
1147	>	//RELATIVE TO the start of the contig
1148
1149	<	/*
1150	<	* remember that text_offset holds the left end of the
1165	<	* alignment, relative to the start of the contig
1166	<	*/
1167	<
1168	<	/*
1169	<	* The 0-based relative position of the left-most character
1170	<	* before the insertion in the contig
1171	<	*/
1149	>	//The 0-based relative position of the left-most character
1150	>	//before the insertion in the contig
1151		int relative_splice_pos = left_splice_pos - atoi(toks[num_extra_toks + left_window_edge_field].c_str());
1152	<	while (pch != NULL)
1153	<	{
1154	<	char* colon = strchr(pch, ':');
1155	<	if (colon)
1156	<	{
1157	<	*colon = 0;
1158	<	int mismatch_pos = atoi(pch);
1159	<
1160	<	/*
1161	<	* for reversely mapped reads,
1162	<	* find the correct mismatched position.
1163	<	*/
1164	<	if (sam_flag & 0x0010){
1165	<	mismatch_pos = spliced_read_len - mismatch_pos - 1;
1166	<	}
1188	<
1189	<	/*
1190	<	* Only count mismatches outside of the insertion region
1191	<	* If there is a mismatch within the insertion,
1192	<	* disallow this hit
1193	<	*/
1194	<	if(mismatch_pos + text_offset <= relative_splice_pos || mismatch_pos + text_offset > relative_splice_pos + insertedSequence.length()){
1195	<	num_mismatches++;
1196	<	}else{
1197	<	return false;
1198	<	}
1152	>	for (int i=0;i<num_mismatches;i++) {
1153	>	int mismatch_pos = mismatches[i];
1154	>	// for reversely mapped reads,
1155	>	//find the correct mismatched position.
1156	>	if (sam_flag & 0x0010){
1157	>	mismatch_pos = spliced_read_len - mismatch_pos - 1;
1158	>	}
1159	>
1160	>	//Only count mismatches outside of the insertion region
1161	>	// If there is a mismatch within the insertion,
1162	>	// disallow this hit
1163	>	if(mismatch_pos + text_offset <= relative_splice_pos || mismatch_pos + text_offset > relative_splice_pos + insertedSequence.length()){
1164	>	spl_num_mismatches++;
1165	>	}else{
1166	>	return false;
1167		}
1200	–	pch = strtok (NULL, ",");
1168		}
1169
1170	<
1171	<	vector<CigarOp> cigar;
1172	<	cigar.push_back(CigarOp(MATCH, left_match_length));
1173	<	cigar.push_back(CigarOp(INS, insertion_match_length));
1174	<	cigar.push_back(CigarOp(MATCH, right_match_length));
1175	<
1170	>	//TODO: merge with the original cigar string
1171	>	vector<CigarOp> splcigar;
1172	>	spliceCigar(splcigar, samcigar, left_match_length, insertion_match_length, right_match_length, INS);
1173	>	/*
1174	>	splcigar.push_back(CigarOp(MATCH, left_match_length));
1175	>	splcigar.push_back(CigarOp(INS, insertion_match_length));
1176	>	splcigar.push_back(CigarOp(MATCH, right_match_length));
1177	>	*/
1178		/*
1179		* For now, disallow hits that don't span
1180		* the insertion. If we allow these types of hits,
#	Line 1219 | Line 1188
1188		bh = create_hit(name,
1189		contig,
1190		left,
1191	<	cigar,
1192	<	sam_flag & 0x0010,  //#######
1191	>	//splcigar,
1192	>	splcigar,
1193	>	sam_flag & 0x0010,
1194		junction_strand == "rev",
1195	<	num_mismatches,
1195	>	spl_num_mismatches,
1196		0,
1197		end);
1198		return true;
1199	<	}
1230	<
1199	>	} //"ins"
1200		else
1201		{
1202		uint32_t left = atoi(toks[num_extra_toks + left_window_edge_field].c_str()) + text_offset;
#	Line 1241 | Line 1210
1210
1211		if ((sam_flag & 0x0010) == 0) //######
1212		{
1213	<	if (left_splice_pos + seg_offset < _anchor_length){
1214	<	return false;
1246	<	}
1213	>	if (left_splice_pos + seg_offset < _anchor_length)
1214	>	return false;
1215		}
1216		else
1217		{
1218	<	if (right_splice_pos + seg_offset < _anchor_length)
1219	<	return false;
1218	>	if (right_splice_pos + seg_offset < _anchor_length)
1219	>	return false;
1220		}
1221		//uint32_t right = atoi(splice_toks[1].c_str()) + right_splice_pos;
1222		//atoi(toks[right_window_edge_field].c_str());
1223		int gap_len = atoi(splice_toks[1].c_str()) - atoi(splice_toks[0].c_str()) - 1;
1224
1225		string junction_strand = toks[num_extra_toks + strand_field];
1226	<	if (!(junction_strand == "rev" || junction_strand == "fwd")||
1227	<	!(orientation == '-' || orientation == '+'))
1226	>	if (!(junction_strand == "rev" || junction_strand == "fwd"))
1227	>	// || !(orientation == '-' || orientation == '+'))
1228		{
1229	<	fprintf(stderr, "Warning: found malformed splice record, skipping\n");
1230	<	//fprintf(stderr, "junction_strand=%s, orientation='%c'\n",
1231	<	//           junction_strand.c_str(), orientation);
1232	<	return false;
1229	>	fprintf(stderr, "Warning: found malformed splice record, skipping\n");
1230	>	//fprintf(stderr, "junction_strand=%s, orientation='%c'\n",
1231	>	//           junction_strand.c_str(), orientation);
1232	>	return false;
1233		}
1234
1267	–	//vector<string> mismatch_toks;
1268	–	char* pch = strtok (mismatches,",");
1235		int mismatches_in_anchor = 0;
1236	<	unsigned char num_mismatches = 0;
1237	<	while (pch != NULL)
1238	<	{
1239	<	char* colon = strchr(pch, ':');
1240	<	if (colon)
1241	<	{
1242	<	*colon = 0;
1243	<	num_mismatches++;
1278	<	int mismatch_pos = atoi(pch);
1279	<	if ((orientation == '+' && abs(mismatch_pos - left_splice_pos) < (int)min_anchor_len) ||
1280	<	(orientation == '-' && abs(((int)spliced_read_len - left_splice_pos + 1) - mismatch_pos)) < (int)min_anchor_len)
1281	<	mismatches_in_anchor++;
1282	<	}
1283	<	//mismatch_toks.push_back(pch);
1284	<	pch = strtok (NULL, ",");
1285	<	}
1236	>	for (int i=0;i<num_mismatches;i++) {
1237	>	spl_num_mismatches++;
1238	>	int mismatch_pos = mismatches[i];
1239	>	if (((sam_flag & 0x0010) == 0 && abs(mismatch_pos - left_splice_pos) < (int)min_anchor_len) ||
1240	>	((sam_flag & 0x0010) != 0 &&
1241	>	abs(((int)spliced_read_len - left_splice_pos + 1) - mismatch_pos)) < (int)min_anchor_len)
1242	>	mismatches_in_anchor++;
1243	>	}
1244
1245		// FIXME: we probably should exclude these hits somewhere, but this
1246		// isn't the right place
1247	<	vector<CigarOp> cigar;
1248	<	cigar.push_back(CigarOp(MATCH, left_splice_pos));
1247	>	vector<CigarOp> splcigar;
1248	>	CigarOpCode opcode=(toks[num_extra_toks + junction_type_field] == "del")? DEL : REF_SKIP ;
1249	>	spliceCigar(splcigar, samcigar, left_splice_pos, gap_len, right_splice_pos, INS);
1250	>	/*
1251	>	splcigar.push_back(CigarOp(MATCH, left_splice_pos));
1252		if(toks[num_extra_toks + junction_type_field] == "del"){
1253	<	cigar.push_back(CigarOp(DEL, gap_len));
1253	>	splcigar.push_back(CigarOp(DEL, gap_len));
1254		}else{
1255	<	cigar.push_back(CigarOp(REF_SKIP, gap_len));
1255	>	splcigar.push_back(CigarOp(REF_SKIP, gap_len));
1256		}
1257	<	cigar.push_back(CigarOp(MATCH, right_splice_pos));
1258	<
1257	>	splcigar.push_back(CigarOp(MATCH, right_splice_pos));
1258	>	*/
1259		bh = create_hit(name,
1260		contig,
1261		left,
1262	<	cigar,
1263	<	orientation == '-',
1262	>	splcigar,
1263	>	(sam_flag & 0x0010),
1264		junction_strand == "rev",
1265	<	num_mismatches,
1265	>	spl_num_mismatches,
1266		mismatches_in_anchor,
1267		end);
1268		return true;
1269		}
1270	<	}
1270	>	} //parse splice data
1271		else
1272		{
1273		fprintf(stderr, "Warning: found malformed splice record, skipping\n");

Diff Legend

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