ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/tophat_cpp/bwt_map.cpp
(Generate patch)
# Line 507 | Line 507
507                          uint32_t left = atoi(toks[num_extra_toks + left_window_edge_field].c_str()) + text_offset;
508                          uint32_t right = left + spliced_read_len - 1;
509  
510 +
511                          /*
512                           * The 0-based position of the last genomic sequence before the insertion
513                           */
# Line 720 | Line 721
721   }
722  
723  
724 < char parseCigar(vector<CigarOp>& cigar, const char* cigar_str,
724 > int parseCigar(vector<CigarOp>& cigar, const char* cigar_str,
725                  bool &spliced_alignment)   {
725 // Mostly pilfered direct from the SAM tools:
726   const char* p_cig = cigar_str;
727 + int refspan=0; //alignment span on reference sequence
728   while (*p_cig)
729   {
730          char* t;
731 <        int length = (int)strtol(p_cig, &t, 10);
732 <        if (length <= 0)
731 >        int op_len = (int)strtol(p_cig, &t, 10);
732 >        if (op_len <= 0)
733          {
734 <                //fprintf (stderr, "CIGAR op has zero length\n");
735 <                return false;
734 >                fprintf (stderr, "Error: CIGAR op has zero length\n");
735 >                return 0;
736          }
737          char op_char = toupper(*t);
738          CigarOpCode opcode;
739 <        if (op_char == 'M') opcode = MATCH;
740 <        else if (op_char == 'I') opcode = INS;
741 <        else if (op_char == 'D') opcode = DEL;
742 <        else if (op_char == 'N')
743 <        {
744 <                if (length > max_report_intron_length)
745 <                        return false;
746 <                opcode = REF_SKIP;
747 <                spliced_alignment = true;
748 <        }
749 <        else if (op_char == 'S') opcode = SOFT_CLIP;
750 <        else if (op_char == 'H') opcode = HARD_CLIP;
751 <        else if (op_char == 'P') opcode = PAD;
752 <        else
753 <        {
754 <                fprintf (stderr, "invalid CIGAR operation\n");
755 <                return false;
756 <        }
739 >        switch (op_char) {
740 >          case '=':
741 >          case 'X':
742 >          case 'M': opcode = MATCH;
743 >                    refspan+=op_len;
744 >                    break;
745 >          case 'I': opcode = INS;
746 >                    break;
747 >    case 'D': opcode = DEL;
748 >              refspan+=op_len;
749 >              break;
750 >    case 'N': if (op_len > max_report_intron_length)
751 >                  return 0;
752 >              opcode = REF_SKIP;
753 >              spliced_alignment = true;
754 >              refspan+=op_len;
755 >              break;
756 >    case 'S': opcode = SOFT_CLIP;
757 >              break;
758 >    case 'H': opcode = HARD_CLIP;
759 >              break;
760 >    case 'P': opcode = PAD;
761 >              break;
762 >    default:  fprintf (stderr, "Error: invalid CIGAR operation\n");
763 >                          return 0;
764 >          }
765          p_cig = t + 1;
766 <        //i += length;
767 <        cigar.push_back(CigarOp(opcode, length));
759 < }
766 >        cigar.push_back(CigarOp(opcode, op_len));
767 > } //while cigar codes
768   if (*p_cig) {
769 <        fprintf (stderr, "unmatched CIGAR operation\n");
770 <        return *p_cig;
769 >          fprintf (stderr, "Error: unmatched CIGAR operation (%s)\n",p_cig);
770 >          return 0;
771      }
772 < return 0;
772 > return refspan;
773   }
774  
775 < void getSAMmismatches(char* &buf, vector<CigarOp>& cigar,
776 <                  int* mismatches, int& num_mismatches, int& sam_nm, bool& antisense_splice) {
777 <        const char* tag_buf = buf;
775 > int getSAMmismatches(char* &buf, vector<CigarOp>& cigar,
776 >                  vector<bool>& mismatches, int& sam_nm, bool& antisense_splice) {
777 >  int gspan=0;//genomic span of the alignment
778 >  const char* tag_buf = buf;
779          sam_nm=0;
780 <        num_mismatches=0;
780 >        int num_mismatches=0;
781          while((tag_buf = get_token((char**)&buf,"\t")))
782          {
783                  vector<string> tuple_fields;
# Line 798 | Line 807
807                    do { p++; } while (isdigit(*p));
808                    bi+=v;
809                    }
810 <                 while (isupper(*p)) {
810 >                 while (isalpha(*p)) {
811                    p++;
803                  mismatches[num_mismatches]=bi;
812                    num_mismatches++;
813 +                  //mismatches.push_back(bi);
814 +                  mismatches[bi]=true;
815                    bi++;
816                    }
817                   if (*p=='^') { //reference deletion
818                     p++;
819 <                   while (isupper(*p)) { //insert read bases
819 >                   while (isalpha(*p)) { //insert read bases
820                            p++; bi++;
821 <                      }
821 >                    }
822                     }
823                  }
824                          }
825 <                        else
826 <                        {
825 >     //else
826 >                        //{
827                                  //fprintf(stderr, "%s attribute not supported\n", tuple_fields[0].c_str());
828                                  //return false;
829 <                        }
829 >                        //}
830                  }
831          }
832           /* By convention,the NM field of the SAM record
# Line 826 | Line 836
836           *  the mismatches due to in/dels
837           */
838          for(vector<CigarOp>::const_iterator itr = cigar.begin(); itr != cigar.end(); ++itr){
839 <                if(itr->opcode == INS){
840 <                        sam_nm -= itr->length;
841 <                }
842 <                if(itr->opcode == DEL){
843 <                        sam_nm -= itr->length;
839 >    switch (itr->opcode)
840 >    {
841 >      case MATCH:
842 >      case REF_SKIP:
843 >      case PAD:
844 >        gspan += itr->length;
845 >        break;
846 >      case DEL:
847 >        gspan += itr->length;
848 >        sam_nm -= itr->length;
849 >        break;
850 >      case INS:
851 >        sam_nm -= itr->length;
852 >        break;
853                  }
854 <        }
854 >         }
855 >        return num_mismatches;
856   }
857  
858   bool SAMHitFactory::get_hit_from_buf(const char* orig_bwt_buf,
# Line 900 | Line 920
920          
921          vector<CigarOp> cigar;
922          bool spliced_alignment = false;
923 <        if (parseCigar(cigar, cigar_str, spliced_alignment)) {
923 >        int refspan=parseCigar(cigar, cigar_str, spliced_alignment);
924 >        if (refspan==0)
925             return false;
905           }
906        
926          //vector<string> attributes;
927          //tokenize(tag_buf, " \t",attributes);
928  
929          bool antisense_splice = false;
911        int num_mismatches = 0;
930          int sam_nm = 0; //the value of the NM tag (edit distance)
931 <        int mismatches[1024];
932 <        getSAMmismatches(buf, cigar, mismatches, num_mismatches, sam_nm, antisense_splice);
931 >        //int mismatches[1024];//array with mismatch positions on the read (0-based from the left aligned end of the read)
932 >        vector<bool> mismatches;
933 >        mismatches.resize(strlen(seq_str), false);
934 >        int num_mismatches=getSAMmismatches(buf, cigar, mismatches, sam_nm, antisense_splice);
935          
936          if (spliced_alignment)
937          {
# Line 941 | Line 961
961          return true;
962   }
963  
964 + void cigar_add(vector<CigarOp>& cigar, CigarOp& op) {
965 + if (op.length<=0) return;
966 + if (cigar.size()>0 && cigar.back().opcode==op.opcode) {
967 +    cigar.back().length+=op.length;
968 +    }
969 + cigar.push_back(op);
970 + }
971  
972 < void spliceCigar(vector<CigarOp>& splcigar, vector<CigarOp>& cigar,
973 <                int8_t l_len, int mlen, int8_t r_len, CigarOpCode opcode) {
972 > int spliceCigar(vector<CigarOp>& splcigar, vector<CigarOp>& cigar, vector<bool> mismatches,
973 >                int &left, int spl_start, int spl_len, CigarOpCode spl_code) {
974   //merge the original 'cigar' with the new insert/gap operation
975 < // at position l_len into splcigar;
976 < //find the original opcode that croses l_len location
977 < int ofs=0;
978 < bool found=false;
979 < for (size_t c = 0 ; c < cigar.size(); ++c)
980 <        {
981 <    ofs+=cigar[c].length;
982 <    if (!found) {
983 <        if (ofs>=l_len) {
984 <           found=true;
985 <           //TODO:inject new code here, splitting existing code if necessary
975 > //at position spl_start and place the result into splcigar;
976 > //TODO: ideally this should also get and rebuild the MD string (alignment mismatches)
977 > int spl_mismatches=0;
978 > //return value: mismatches in the insert region for INS case,
979 > //or number of mismatches in the anchor region
980 > //return -1 if somehow the hit seems bad
981 >
982 > //these offsets are relative to the beginning of alignment on reference
983 > int spl_ofs=spl_start-left; //relative position of splice op
984 > int spl_ofs_end=spl_ofs; //relative position of first ref base AFTER splice op
985 > CigarOp gapop(spl_code, spl_len); //for DEL or REF_SKIP
986 > if (spl_code==INS)
987 >     spl_ofs_end += spl_len;
988 > int ref_ofs=0; //working offset on reference
989 > int read_ofs=0; //working offset on the read, relative to the leftmost aligned base
990 > bool xfound=false;
991 > //if (left<=spl_start+spl_len) {
992 > if (spl_ofs_end>0) {
993 >     int prev_opcode=0;
994 >     int prev_oplen=0;
995 >     for (size_t c = 0 ; c < cigar.size(); ++c)
996 >      {
997 >        int prev_read_ofs=read_ofs;
998 >        int cur_op_ofs=ref_ofs;
999 >        int cur_opcode=cigar[c].opcode;
1000 >        int cur_oplen=cigar[c].length;
1001 >        switch (cur_opcode) {
1002 >           case MATCH:
1003 >             ref_ofs+=cur_oplen;
1004 >             read_ofs+=cur_oplen;
1005 >             if (spl_code==REF_SKIP)
1006 >                for (int o=cur_op_ofs;o<ref_ofs;o++) {
1007 >                    int rofs=prev_read_ofs+(cur_op_ofs-o);
1008 >                    if (abs(spl_ofs-o)<min_anchor_len && mismatches[rofs])
1009 >                      spl_mismatches++;
1010 >                    }
1011 >             break;
1012 >           case DEL:
1013 >           case REF_SKIP:
1014 >           case PAD:
1015 >             ref_ofs+=cur_oplen;
1016 >             break;
1017 >           case SOFT_CLIP:
1018 >           case INS:
1019 >             read_ofs+=cur_oplen;
1020 >             break;
1021 >           //case HARD_CLIP:
1022 >           }
1023  
1024 <           continue;
1024 >        if (cur_op_ofs>=spl_ofs_end || ref_ofs<=spl_ofs) {
1025 >           if (cur_op_ofs==spl_ofs_end && spl_code!=INS)
1026 >              {
1027 >              xfound=true;
1028 >               //we have to insert the gap here first
1029 >              cigar_add(splcigar, gapop);
1030 >              //also, check
1031 >              }
1032 >           cigar_add(splcigar, cigar[c]);
1033             }
1034 <        else {
1035 <           //not found yet, just copy these codes
1036 <           splcigar.push_back(cigar[c]);
1037 <           }
1038 <       }
1039 <    else { //found already
1040 <       //TODO: check this
1041 <       splcigar.push_back(cigar[c]);
1042 <       }
1043 <        }
1034 >        else //if (ref_ofs>spl_ofs) {
1035 >           { //op intersection
1036 >           xfound=true;
1037 >           if (spl_code==INS) {
1038 >                 //we have to shorten cur_opcode
1039 >                 // find the overlap between current range
1040 >                 int ovl_start = (cur_op_ofs>spl_ofs) ? cur_op_ofs : spl_ofs;
1041 >                 int ovl_end = (ref_ofs>spl_ofs_end) ? spl_ofs_end : ref_ofs;
1042 >                 CigarOp op(cigar[c]);
1043 >                 int len=op.length;
1044 >                 len-=(ovl_end-ovl_start);
1045 >                 if (len>0) {
1046 >                     op.length=len;
1047 >                     cigar_add(splcigar, op);
1048 >                     }
1049 >                 }
1050 >              else {//DEL or REF_SKIP
1051 >                 //spl_ofs == spl_ofs_end
1052 >                 //we have to split cur_opcode
1053 >                 //look for mismatches within min_anchor_len distance from splice point
1054 >                 CigarOp op(cigar[c]);
1055 >                 op.length=spl_ofs-cur_op_ofs;
1056 >                 cigar_add(splcigar, op);
1057 >                 cigar_add(splcigar, gapop);
1058 >                 op.length=ref_ofs-spl_ofs;
1059 >                 cigar_add(splcigar,op);
1060 >                 }
1061 >            } //op intersection
1062 >        prev_opcode=cur_opcode;
1063 >        prev_oplen=cur_oplen;
1064 >      } //for each cigar opcode
1065 >    } //intersection possible
1066 >
1067 >  //if (!xfound) {//no intersection found between splice event and alignment
1068 >   if (spl_ofs_end<=0) {
1069 >      //alignment starts after the splice event
1070 >      if (spl_code==INS) left-=spl_len;
1071 >                   else  left+=spl_len;
1072 >      }
1073 >   //else {
1074 >     //alignment ends before the splice event
1075 >     //nothing to do
1076 >    //  }
1077 >   //return spl_mismatches;
1078 >  // }
1079 >
1080 > return spl_mismatches;
1081   }
1082  
1083   bool SplicedSAMHitFactory::get_hit_from_buf(const char* orig_bwt_buf,
# Line 998 | Line 1107
1107          const char* mate_ref_str =  get_token((char**)&buf,"\t");
1108          const char* mate_pos_str =  get_token((char**)&buf,"\t");
1109          const char* inferred_insert_sz_str =  get_token((char**)&buf,"\t");
1110 <        int num_mismatches=0;
1111 <    int mismatches[1024]; //list of 0-based mismatch positions in this read
1110 >        //int num_mismatches=0;
1111 >  //int mismatches[1024]; //list of 0-based mismatch positions in this read
1112                            //parsed from SAM's MD:Z: tag
1113          const char* seq_str =  get_token((char**)&buf,"\t");
1114          if (seq)
# Line 1028 | Line 1137
1137  
1138          int sam_flag = atoi(sam_flag_str);
1139          int text_offset = atoi(text_offset_str);
1140 <    text_offset--; //SAM is 1-based, Bowtie is 0-based
1140 >  text_offset--; //make it 0-based (SAM is 1-based, Bowtie is 0-based)
1141          bool end = true;
1142          unsigned int seg_offset = 0;
1143          unsigned int seg_num = 0;
# Line 1039 | Line 1148
1148  
1149          vector<CigarOp> samcigar;
1150          bool spliced_alignment = false;
1151 <        if (parseCigar(samcigar, cigar_str, spliced_alignment)) {
1152 <           return false;
1153 <           }
1045 <
1151 >        int refspan=parseCigar(samcigar, cigar_str, spliced_alignment);
1152 >  if (refspan==0)
1153 >      return false;
1154          bool antisense_splice = false;
1155          int sam_nm = 0;
1156 <    getSAMmismatches(buf, samcigar, mismatches, num_mismatches, sam_nm, antisense_splice);
1156 >  vector<bool> mismatches;
1157 >  mismatches.resize(strlen(seq_str), false);
1158 >
1159 >  int num_mismatches=getSAMmismatches(buf, samcigar, mismatches, sam_nm, antisense_splice);
1160  
1161          //##############################################
1162  
# Line 1054 | Line 1165
1165  
1166          // Parse the text_name field to recover the splice coords
1167          vector<string> toks;
1057
1168          tokenize_strict(text_name, "|", toks);
1169  
1170          int num_extra_toks = (int)toks.size() - 6;
# Line 1083 | Line 1193
1193                          //        toks[num_extra_toks + splice_field].c_str());
1194                          return false;
1195                  }
1196 <        int spl_num_mismatches=0;
1196 >
1197 >    string junction_strand = toks[num_extra_toks + strand_field];
1198 >    if(junction_strand != "rev" && junction_strand != "fwd"){
1199 >      fprintf(stderr, "Malformed insertion record\n");
1200 >      return false;
1201 >      }
1202 >
1203                  //
1204                  // check for an insertion hit
1205                  //
1206                  if(toks[num_extra_toks + junction_type_field] == "ins")
1207                    {
1208 <                        int8_t spliced_read_len = strlen(seq_str);
1208 >                        //int8_t spliced_read_len = strlen(seq_str);
1209 >                        //TODO FIXME: use the CIGAR instead of seq length!
1210                           // The 0-based position of the left edge of the alignment. Note that this
1211 <                         // value may need to be futher corrected to account for the presence of
1211 >                 // value may need to be further corrected to account for the presence of
1212                           // of the insertion.
1213 <                        uint32_t left = atoi(toks[num_extra_toks + left_window_edge_field].c_str()) + text_offset;
1097 <                        uint32_t right = left + spliced_read_len - 1;
1213 >                        int left = atoi(toks[num_extra_toks + left_window_edge_field].c_str()) + text_offset;
1214  
1215                          // The 0-based position of the last genomic sequence before the insertion
1216 <                        uint32_t left_splice_pos = atoi(splice_toks[0].c_str());
1216 >                        int left_splice_pos = atoi(splice_toks[0].c_str());
1217  
1218                          string insertedSequence = splice_toks[1];
1219                          // The 0-based position of the first genomic sequence after the insertion
1220  
1221 +      vector<CigarOp> splcigar;
1222 +      //this also updates left to the adjusted genomic coordinates
1223 +      int spl_num_mismatches=spliceCigar(splcigar, samcigar, mismatches,
1224 +                        left, left_splice_pos+1, insertedSequence.length(), INS);
1225 +      if (spl_num_mismatches<0) return false;
1226 +      num_mismatches-=spl_num_mismatches;
1227 +      /*
1228                          uint32_t right_splice_pos = left_splice_pos + 1;
1229 +
1230 +      //uint32_t right = left + spliced_read_len - 1;
1231 +      int right = left + refspan - 1;
1232 +
1233                          if(left > left_splice_pos){
1234 <                                /*
1235 <                                 * The genomic position of the left edge of the alignment needs to be corrected
1236 <                                 * If the alignment does not extend into the insertion, simply subtract the length
1110 <                                 * of the inserted sequence, otherwise, just set it equal to the right edge
1111 <                                 */
1234 >                                 //The genomic position of the left edge of the alignment needs to be corrected
1235 >                                 //If the alignment does not extend into the insertion, simply subtract the length
1236 >                                 //of the inserted sequence, otherwise, just set it equal to the right edge
1237                                  left = left - insertedSequence.length();
1238                                  if(left < right_splice_pos){
1239                                          left = right_splice_pos;
# Line 1120 | Line 1245
1245                                          right = left_splice_pos;
1246                                  }
1247                          }
1248 <                        /*
1249 <                         * Now, right and left should be properly transformed into genomic coordinates
1250 <                         * We should be able to deduce how much the alignment matches the insertion
1126 <                         * simply based on the length of the read
1127 <                         */
1248 >                 // Now, right and left should be properly transformed into genomic coordinates
1249 >                 // We should be able to deduce how much the alignment matches the insertion
1250 >                 // simply based on the length of the read
1251                          int left_match_length = 0;
1252                          if(left <= left_splice_pos){
1253                                  left_match_length = left_splice_pos - left + 1;
# Line 1138 | Line 1261
1261                          if(left_match_length <= 0 || right_match_length <= 0 || insertion_match_length <= 0)
1262                            return false;
1263  
1141                        string junction_strand = toks[num_extra_toks + strand_field];
1142                        if(junction_strand != "rev" && junction_strand != "fwd"){
1143                                fprintf(stderr, "Malformed insertion record\n");
1144                                return false;
1145                        }
1146
1264                          //char* pch = strtok( mismatches, ",");
1265                          //unsigned char num_mismatches = 0;
1266                          //text_offset holds the left end of the alignment,
# Line 1152 | Line 1269
1269                          //The 0-based relative position of the left-most character
1270                          //before the insertion in the contig
1271                          int relative_splice_pos = left_splice_pos - atoi(toks[num_extra_toks + left_window_edge_field].c_str());
1272 <                        for (int i=0;i<num_mismatches;i++) {
1272 >                        for (size_t i=0;i<mismatches.size();++i) {
1273                                  int mismatch_pos = mismatches[i];
1274                                  // for reversely mapped reads,
1275                                  //find the correct mismatched position.
# Line 1163 | Line 1280
1280                                  //Only count mismatches outside of the insertion region
1281                                  // If there is a mismatch within the insertion,
1282                                  // disallow this hit
1283 <                                if(mismatch_pos + text_offset <= relative_splice_pos || mismatch_pos + text_offset > relative_splice_pos + insertedSequence.length()){
1283 >                                if(mismatch_pos + text_offset <= relative_splice_pos ||
1284 >                                    mismatch_pos + text_offset > relative_splice_pos + insertedSequence.length()){
1285                                    spl_num_mismatches++;
1286                                  }else{
1287                                    return false;
1288                                  }
1289                          }
1290 +      */
1291 +                        //vector<CigarOp> splcigar;
1292 +                        //spliceCigar(splcigar, samcigar, left_match_length, insertion_match_length, right_match_length, INS);
1293 +                        //splcigar.push_back(CigarOp(MATCH, left_match_length));
1294 +                        //splcigar.push_back(CigarOp(INS, insertion_match_length));
1295 +                        //splcigar.push_back(CigarOp(MATCH, right_match_length));
1296  
1173            //TODO: merge with the original cigar string
1174                        vector<CigarOp> splcigar;
1175                        spliceCigar(splcigar, samcigar, left_match_length, insertion_match_length, right_match_length, INS);
1176                        /*
1177                        splcigar.push_back(CigarOp(MATCH, left_match_length));
1178                        splcigar.push_back(CigarOp(INS, insertion_match_length));
1179                        splcigar.push_back(CigarOp(MATCH, right_match_length));
1180            */
1181                        /*
1182                         * For now, disallow hits that don't span
1183                         * the insertion. If we allow these types of hits,
1184                         * then long_spanning.cpp needs to be updated
1185                         * in order to intelligently merge these kinds
1186                         * of reads back together
1187                         *
1188                         * Following code has been changed to allow segment that end
1189                         * in an insertion
1190                         */
1297                          bh = create_hit(name,
1298                                          contig,
1299                                          left,
# Line 1195 | Line 1301
1301                                          splcigar,
1302                                          sam_flag & 0x0010,
1303                                          junction_strand == "rev",
1304 <                                        spl_num_mismatches,
1304 >                                        num_mismatches,
1305                                          0,
1306                                          end);
1307                          return true;
1308                  } //"ins"
1309 <                else
1309 >                else //"del" or intron
1310                    {
1311 <                    uint32_t left = atoi(toks[num_extra_toks + left_window_edge_field].c_str()) + text_offset;
1312 <                    int spliced_read_len = strlen(seq_str);
1207 <                    int8_t left_splice_pos = atoi(splice_toks[0].c_str()) - left + 1;
1208 <                    if(left_splice_pos > spliced_read_len) left_splice_pos = spliced_read_len;
1209 <                    int8_t right_splice_pos = spliced_read_len - left_splice_pos;
1311 >                    // The 0-based position of the left edge of the alignment.
1312 >                    int left = atoi(toks[num_extra_toks + left_window_edge_field].c_str()) + text_offset;
1313  
1314 <                    if (right_splice_pos <= 0 || left_splice_pos <= 0)
1315 <                      return false;
1314 >                    // The 0-based position of the last genomic sequence before the deletion
1315 >              int left_splice_pos = atoi(splice_toks[0].c_str());
1316  
1317 +        int gap_len = atoi(splice_toks[1].c_str()) - left_splice_pos - 1;
1318 +        /*
1319                      if ((sam_flag & 0x0010) == 0) //######
1320                        {
1321 <                                if (left_splice_pos + seg_offset < _anchor_length)
1322 <                                  return false;
1323 <                      }
1324 <                    else
1325 <                      {
1326 <                                if (right_splice_pos + seg_offset < _anchor_length)
1327 <                                  return false;
1328 <                      }
1321 >            if (left_splice_ofs + seg_offset < _anchor_length)
1322 >              return false;
1323 >          }
1324 >        else
1325 >          {
1326 >            if (right_splice_ofs + seg_offset < _anchor_length)
1327 >              return false;
1328 >          }
1329 >          */
1330                      //uint32_t right = atoi(splice_toks[1].c_str()) + right_splice_pos;
1331                      //atoi(toks[right_window_edge_field].c_str());
1226                    int gap_len = atoi(splice_toks[1].c_str()) - atoi(splice_toks[0].c_str()) - 1;
1332  
1333 <                    string junction_strand = toks[num_extra_toks + strand_field];
1334 <                    if (!(junction_strand == "rev" || junction_strand == "fwd"))
1335 <                        // || !(orientation == '-' || orientation == '+'))
1231 <                      {
1232 <                                fprintf(stderr, "Warning: found malformed splice record, skipping\n");
1233 <                                //fprintf(stderr, "junction_strand=%s, orientation='%c'\n",
1234 <                                //           junction_strand.c_str(), orientation);
1235 <                                return false;
1236 <                      }
1333 >        /*
1334 >        //offset of deletion point, relative to the beginning of the alignment
1335 >        int left_splice_ofs = left_splice_pos-left+1;
1336  
1337                      int mismatches_in_anchor = 0;
1338 <                    for (int i=0;i<num_mismatches;i++) {
1339 <                                spl_num_mismatches++;
1340 <                                int mismatch_pos = mismatches[i];
1341 <                                if (((sam_flag & 0x0010) == 0 && abs(mismatch_pos - left_splice_pos) < (int)min_anchor_len) ||
1342 <                                  ((sam_flag & 0x0010) != 0 &&
1343 <                                                  abs(((int)spliced_read_len - left_splice_pos + 1) - mismatch_pos)) < (int)min_anchor_len)
1338 >                    for (size_t i=0;i<mismatches.size();++i) {
1339 >                                  spl_num_mismatches++;
1340 >                                  int mismatch_pos = mismatches[i];
1341 >                                  if (((sam_flag & 0x0010) == 0 && abs(mismatch_pos - left_splice_ofs) < (int)min_anchor_len) ||
1342 >                                    ((sam_flag & 0x0010) != 0 &&
1343 >                                                  abs(((int)refspan - left_splice_ofs + 1) - mismatch_pos)) < (int)min_anchor_len)
1344                                            mismatches_in_anchor++;
1345 <                                 }
1346 <
1248 <                    // FIXME: we probably should exclude these hits somewhere, but this
1249 <                    // isn't the right place
1345 >                                  }
1346 >         */
1347                      vector<CigarOp> splcigar;
1348                      CigarOpCode opcode=(toks[num_extra_toks + junction_type_field] == "del")? DEL : REF_SKIP ;
1349 <                        spliceCigar(splcigar, samcigar, left_splice_pos, gap_len, right_splice_pos, INS);
1350 <            /*
1349 >
1350 >                    int spl_num_mismatches=spliceCigar(splcigar, samcigar, mismatches, left,
1351 >                        left_splice_pos+1, gap_len, opcode);
1352 >                    if (spl_num_mismatches<0) // || spl_num_mismatches>max_anchor_mismatches)
1353 >                         return false;
1354 >        /*
1355                      splcigar.push_back(CigarOp(MATCH, left_splice_pos));
1356                      if(toks[num_extra_toks + junction_type_field] == "del"){
1357                        splcigar.push_back(CigarOp(DEL, gap_len));
# Line 1258 | Line 1359
1359                        splcigar.push_back(CigarOp(REF_SKIP, gap_len));
1360                      }
1361                      splcigar.push_back(CigarOp(MATCH, right_splice_pos));
1362 <            */
1362 >        */
1363                      bh = create_hit(name,
1364                                      contig,
1365                                      left,
1366                                      splcigar,
1367                                      (sam_flag & 0x0010),
1368                                      junction_strand == "rev",
1369 +                                    num_mismatches,
1370                                      spl_num_mismatches,
1269                                    mismatches_in_anchor,
1371                                      end);
1372                      return true;
1373                    }
# Line 1283 | Line 1384
1384   }
1385  
1386  
1387 +
1388   bool BAMHitFactory::get_hit_from_buf(const char* orig_bwt_buf,
1389                                   BowtieHit& bh, bool strip_slash,
1390                                  char* name_out, char* name_tags,
# Line 1600 | Line 1702
1702                          case MATCH:
1703                                  for (size_t m = 0; m < cigar[c].length; ++m)
1704                                  {
1705 <                                        if (DoC[j + m] < 0xFFFFFFFF)
1705 >                                        if (DoC[j + m] < VMAXINT32)
1706                                                  DoC[j + m]++;
1707                                  }
1708                                  //fall through this case to REF_SKIP is intentional

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines