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<int>& 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 >        //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++;
812 <                  mismatches[num_mismatches]=bi;
804 <                  num_mismatches++;
812 >                  mismatches.push_back(bi);
813                    bi++;
814                    }
815                   if (*p=='^') { //reference deletion
816                     p++;
817 <                   while (isupper(*p)) { //insert read bases
817 >                   while (isalpha(*p)) { //insert read bases
818                            p++; bi++;
819 <                      }
819 >                    }
820                     }
821                  }
822                          }
# Line 826 | Line 834
834           *  the mismatches due to in/dels
835           */
836          for(vector<CigarOp>::const_iterator itr = cigar.begin(); itr != cigar.end(); ++itr){
837 <                if(itr->opcode == INS){
838 <                        sam_nm -= itr->length;
839 <                }
840 <                if(itr->opcode == DEL){
841 <                        sam_nm -= itr->length;
837 >    switch (itr->opcode)
838 >    {
839 >      case MATCH:
840 >      case REF_SKIP:
841 >      case PAD:
842 >        gspan += itr->length;
843 >        break;
844 >      case DEL:
845 >        gspan += itr->length;
846 >        sam_nm -= itr->length;
847 >        break;
848 >      case INS:
849 >        sam_nm -= itr->length;
850 >        break;
851                  }
852 <        }
852 >         }
853 >        return gspan;
854   }
855  
856   bool SAMHitFactory::get_hit_from_buf(const char* orig_bwt_buf,
# Line 900 | Line 918
918          
919          vector<CigarOp> cigar;
920          bool spliced_alignment = false;
921 <        if (parseCigar(cigar, cigar_str, spliced_alignment)) {
921 >        int refspan=parseCigar(cigar, cigar_str, spliced_alignment);
922 >        if (refspan==0)
923             return false;
905           }
906        
924          //vector<string> attributes;
925          //tokenize(tag_buf, " \t",attributes);
926  
927          bool antisense_splice = false;
911        int num_mismatches = 0;
928          int sam_nm = 0; //the value of the NM tag (edit distance)
929 <        int mismatches[1024];
930 <        getSAMmismatches(buf, cigar, mismatches, num_mismatches, sam_nm, antisense_splice);
929 >        //int mismatches[1024];//array with mismatch positions on the read (0-based from the left aligned end of the read)
930 >        vector<int> mismatches;
931 >        getSAMmismatches(buf, cigar, mismatches, sam_nm, antisense_splice);
932          
933          if (spliced_alignment)
934          {
# Line 921 | Line 938
938                                  cigar,
939                                  sam_flag & 0x0010,
940                                  antisense_splice,
941 <                                num_mismatches,
941 >                                mismatches.size(),
942                                  0,
943                                  end);
944          }
# Line 934 | Line 951
951                                  cigar,
952                                  sam_flag & 0x0010,
953                                  false,
954 <                                num_mismatches,
954 >                                mismatches.size(),
955                                  0,
956                                  end);
957          }
958          return true;
959   }
960  
961 + void cigar_add(vector<CigarOp>& cigar, CigarOp& op) {
962 + if (op.length<=0) return;
963 + if (cigar.size()>0 && cigar.back().opcode==op.opcode) {
964 +    cigar.back().length+=op.length;
965 +    }
966 + cigar.push_back(op);
967 + }
968  
969 < void spliceCigar(vector<CigarOp>& splcigar, vector<CigarOp>& cigar,
970 <                int8_t l_len, int mlen, int8_t r_len, CigarOpCode opcode) {
969 > int spliceCigar(vector<CigarOp>& splcigar, vector<CigarOp>& cigar, vector<int> mismatches,
970 >                int &left, int spl_start, int spl_len, CigarOpCode spl_code) {
971   //merge the original 'cigar' with the new insert/gap operation
972 < // at position l_len into splcigar;
973 < //find the original opcode that croses l_len location
974 < int ofs=0;
975 < bool found=false;
976 < for (size_t c = 0 ; c < cigar.size(); ++c)
977 <        {
978 <    ofs+=cigar[c].length;
979 <    if (!found) {
980 <        if (ofs>=l_len) {
981 <           found=true;
982 <           //TODO:inject new code here, splitting existing code if necessary
983 <
984 <           continue;
972 > //at position spl_start and place the result into splcigar;
973 > int num_mismatches=mismatches.size();
974 > //these offsets are relative to the beginning of alignment on reference
975 > int spl_ofs=spl_start-left; //relative position of splice op
976 > int spl_ofs_end=spl_ofs; //relative position of first ref base AFTER splice op
977 > CigarOp gapop(spl_code, spl_len); //for DEL or REF_SKIP
978 > if (spl_code==INS)
979 >     spl_ofs_end += spl_len;
980 > int ref_ofs=0; //working offset on reference
981 > int read_ofs=0; //working offset on the read, relative to the leftmost aligned base
982 > bool xfound=false;
983 > //if (left<=spl_start+spl_len) {
984 > if (spl_ofs_end>0) {
985 >     int prev_opcode=0;
986 >     int prev_oplen=0;
987 >     for (size_t c = 0 ; c < cigar.size(); ++c)
988 >      {
989 >        int prev_read_ofs=read_ofs;
990 >        int cur_op_ofs=ref_ofs;
991 >        int cur_opcode=cigar[c].opcode;
992 >        int cur_oplen=cigar[c].length;
993 >        switch (cur_opcode) {
994 >           case MATCH:
995 >             ref_ofs+=cur_oplen;
996 >             read_ofs+=cur_oplen;
997 >             break;
998 >           case DEL:
999 >           case REF_SKIP:
1000 >           case PAD:
1001 >             ref_ofs+=cur_oplen;
1002 >             break;
1003 >           case SOFT_CLIP:
1004 >           case INS:
1005 >             read_ofs+=cur_oplen;
1006 >             break;
1007 >           //case HARD_CLIP:
1008             }
1009 <        else {
1010 <           //not found yet, just copy these codes
1011 <           splcigar.push_back(cigar[c]);
1012 <           }
1013 <       }
1014 <    else { //found already
1015 <       //TODO: check this
1016 <       splcigar.push_back(cigar[c]);
1017 <       }
1018 <        }
1009 >        if (cur_op_ofs>=spl_ofs_end || ref_ofs<=spl_ofs) {
1010 >           if (spl_code!=INS && cur_op_ofs==spl_ofs_end)
1011 >              {
1012 >               //we have to insert the gap here first
1013 >               cigar_add(splcigar, gapop);
1014 >              }
1015 >           cigar_add(splcigar, cigar[c]);
1016 >           }
1017 >        else //if (ref_ofs>spl_ofs) {
1018 >           { //op intersection
1019 >           xfound=true;
1020 >           if (spl_code==INS) {
1021 >                 //we have to shorten cur_opcode
1022 >                 // find the overlap between current range
1023 >                 int ovl_start = (cur_op_ofs>spl_ofs) ? cur_op_ofs : spl_ofs;
1024 >                 int ovl_end = (ref_ofs>spl_ofs_end) ? spl_ofs_end : ref_ofs;
1025 >                 CigarOp op(cigar[c]);
1026 >                 int len=op.length;
1027 >                 len-=(ovl_end-ovl_start);
1028 >                 if (len>0) {
1029 >                     op.length=len;
1030 >                     cigar_add(splcigar, op);
1031 >                     }
1032 >                 }
1033 >              else {//DEL or REF_SKIP
1034 >                 //spl_ofs == spl_ofs_end
1035 >                 //we have to split cur_opcode
1036 >                 CigarOp op(cigar[c]);
1037 >                 op.length=spl_ofs-cur_op_ofs;
1038 >                 cigar_add(splcigar, op);
1039 >                 cigar_add(splcigar, gapop);
1040 >                 op.length=ref_ofs-spl_ofs;
1041 >                 cigar_add(splcigar,op);
1042 >                 }
1043 >            } //op intersection
1044 >        prev_opcode=cur_opcode;
1045 >        prev_oplen=cur_oplen;
1046 >      } //for each cigar opcode
1047 >    } //intersection possible
1048 >  if (!xfound) {//no intersection found between splice event and alignment
1049 >   if (spl_ofs_end<=0) {
1050 >      //alignment starts after the splice event
1051 >      if (spl_code==INS) left-=spl_len;
1052 >                   else  left+=spl_len;
1053 >      }
1054 >   //else {
1055 >     //alignment ends before the splice event
1056 >     //nothing to do
1057 >    //  }
1058 >   }
1059 >
1060 > //TODO: adjust mismatches if opcode==INS and they fall in the removed part
1061 > //TODO: what about anchor mismatches for REF_SKIP
1062 > return num_mismatches;
1063   }
1064  
1065   bool SplicedSAMHitFactory::get_hit_from_buf(const char* orig_bwt_buf,
# Line 998 | Line 1089
1089          const char* mate_ref_str =  get_token((char**)&buf,"\t");
1090          const char* mate_pos_str =  get_token((char**)&buf,"\t");
1091          const char* inferred_insert_sz_str =  get_token((char**)&buf,"\t");
1092 <        int num_mismatches=0;
1093 <    int mismatches[1024]; //list of 0-based mismatch positions in this read
1092 >        //int num_mismatches=0;
1093 >  //int mismatches[1024]; //list of 0-based mismatch positions in this read
1094                            //parsed from SAM's MD:Z: tag
1095          const char* seq_str =  get_token((char**)&buf,"\t");
1096          if (seq)
# Line 1028 | Line 1119
1119  
1120          int sam_flag = atoi(sam_flag_str);
1121          int text_offset = atoi(text_offset_str);
1122 <    text_offset--; //SAM is 1-based, Bowtie is 0-based
1122 >  text_offset--; //make it 0-based (SAM is 1-based, Bowtie is 0-based)
1123          bool end = true;
1124          unsigned int seg_offset = 0;
1125          unsigned int seg_num = 0;
# Line 1039 | Line 1130
1130  
1131          vector<CigarOp> samcigar;
1132          bool spliced_alignment = false;
1133 <        if (parseCigar(samcigar, cigar_str, spliced_alignment)) {
1134 <           return false;
1135 <           }
1045 <
1133 >        int refspan=parseCigar(samcigar, cigar_str, spliced_alignment);
1134 >  if (refspan==0)
1135 >      return false;
1136          bool antisense_splice = false;
1137          int sam_nm = 0;
1138 <  getSAMmismatches(buf, samcigar, mismatches, num_mismatches, sam_nm, antisense_splice);
1138 >  vector<int> mismatches;
1139 >  getSAMmismatches(buf, samcigar, mismatches, sam_nm, antisense_splice);
1140  
1141          //##############################################
1142  
# Line 1054 | Line 1145
1145  
1146          // Parse the text_name field to recover the splice coords
1147          vector<string> toks;
1057
1148          tokenize_strict(text_name, "|", toks);
1149  
1150          int num_extra_toks = (int)toks.size() - 6;
# Line 1083 | Line 1173
1173                          //        toks[num_extra_toks + splice_field].c_str());
1174                          return false;
1175                  }
1176 <        int spl_num_mismatches=0;
1176 >
1177 >    string junction_strand = toks[num_extra_toks + strand_field];
1178 >    if(junction_strand != "rev" && junction_strand != "fwd"){
1179 >      fprintf(stderr, "Malformed insertion record\n");
1180 >      return false;
1181 >      }
1182 >
1183 >    int spl_num_mismatches=0;
1184 >
1185                  //
1186                  // check for an insertion hit
1187                  //
1188                  if(toks[num_extra_toks + junction_type_field] == "ins")
1189                    {
1190 <                        int8_t spliced_read_len = strlen(seq_str);
1190 >                        //int8_t spliced_read_len = strlen(seq_str);
1191 >                        //TODO FIXME: use the CIGAR instead of seq length!
1192                           // The 0-based position of the left edge of the alignment. Note that this
1193 <                         // value may need to be futher corrected to account for the presence of
1193 >                 // value may need to be further corrected to account for the presence of
1194                           // of the insertion.
1195 <                        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;
1195 >                        int left = atoi(toks[num_extra_toks + left_window_edge_field].c_str()) + text_offset;
1196  
1197                          // The 0-based position of the last genomic sequence before the insertion
1198 <                        uint32_t left_splice_pos = atoi(splice_toks[0].c_str());
1198 >                        int left_splice_pos = atoi(splice_toks[0].c_str());
1199  
1200                          string insertedSequence = splice_toks[1];
1201                          // The 0-based position of the first genomic sequence after the insertion
1202  
1203 +      vector<CigarOp> splcigar;
1204 +
1205 +      //TODO: implement this
1206 +      //this also updates left to the adjusted genomic coordinates
1207 +      spl_num_mismatches=spliceCigar(splcigar, samcigar, mismatches,
1208 +                        left, left_splice_pos, insertedSequence.length(), INS);
1209 +      /*
1210                          uint32_t right_splice_pos = left_splice_pos + 1;
1211 +
1212 +      //uint32_t right = left + spliced_read_len - 1;
1213 +      int right = left + refspan - 1;
1214 +
1215                          if(left > left_splice_pos){
1216 <                                /*
1217 <                                 * The genomic position of the left edge of the alignment needs to be corrected
1218 <                                 * 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 <                                 */
1216 >                                 //The genomic position of the left edge of the alignment needs to be corrected
1217 >                                 //If the alignment does not extend into the insertion, simply subtract the length
1218 >                                 //of the inserted sequence, otherwise, just set it equal to the right edge
1219                                  left = left - insertedSequence.length();
1220                                  if(left < right_splice_pos){
1221                                          left = right_splice_pos;
# Line 1120 | Line 1227
1227                                          right = left_splice_pos;
1228                                  }
1229                          }
1230 <                        /*
1231 <                         * Now, right and left should be properly transformed into genomic coordinates
1232 <                         * We should be able to deduce how much the alignment matches the insertion
1126 <                         * simply based on the length of the read
1127 <                         */
1230 >                 // Now, right and left should be properly transformed into genomic coordinates
1231 >                 // We should be able to deduce how much the alignment matches the insertion
1232 >                 // simply based on the length of the read
1233                          int left_match_length = 0;
1234                          if(left <= left_splice_pos){
1235                                  left_match_length = left_splice_pos - left + 1;
# Line 1138 | Line 1243
1243                          if(left_match_length <= 0 || right_match_length <= 0 || insertion_match_length <= 0)
1244                            return false;
1245  
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
1246                          //char* pch = strtok( mismatches, ",");
1247                          //unsigned char num_mismatches = 0;
1248                          //text_offset holds the left end of the alignment,
# Line 1152 | Line 1251
1251                          //The 0-based relative position of the left-most character
1252                          //before the insertion in the contig
1253                          int relative_splice_pos = left_splice_pos - atoi(toks[num_extra_toks + left_window_edge_field].c_str());
1254 <                        for (int i=0;i<num_mismatches;i++) {
1254 >                        for (size_t i=0;i<mismatches.size();++i) {
1255                                  int mismatch_pos = mismatches[i];
1256                                  // for reversely mapped reads,
1257                                  //find the correct mismatched position.
# Line 1163 | Line 1262
1262                                  //Only count mismatches outside of the insertion region
1263                                  // If there is a mismatch within the insertion,
1264                                  // disallow this hit
1265 <                                if(mismatch_pos + text_offset <= relative_splice_pos || mismatch_pos + text_offset > relative_splice_pos + insertedSequence.length()){
1265 >                                if(mismatch_pos + text_offset <= relative_splice_pos ||
1266 >                                    mismatch_pos + text_offset > relative_splice_pos + insertedSequence.length()){
1267                                    spl_num_mismatches++;
1268                                  }else{
1269                                    return false;
1270                                  }
1271                          }
1272 +      */
1273 +                        //vector<CigarOp> splcigar;
1274 +                        //spliceCigar(splcigar, samcigar, left_match_length, insertion_match_length, right_match_length, INS);
1275 +                        //splcigar.push_back(CigarOp(MATCH, left_match_length));
1276 +                        //splcigar.push_back(CigarOp(INS, insertion_match_length));
1277 +                        //splcigar.push_back(CigarOp(MATCH, right_match_length));
1278  
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                         */
1279                          bh = create_hit(name,
1280                                          contig,
1281                                          left,
# Line 1200 | Line 1288
1288                                          end);
1289                          return true;
1290                  } //"ins"
1291 <                else
1291 >                else //"del" or intron
1292                    {
1293 <                    uint32_t left = atoi(toks[num_extra_toks + left_window_edge_field].c_str()) + text_offset;
1294 <                    int spliced_read_len = strlen(seq_str);
1295 <                    int8_t left_splice_pos = atoi(splice_toks[0].c_str()) - left + 1;
1293 >                    int left = atoi(toks[num_extra_toks + left_window_edge_field].c_str()) + text_offset;
1294 >                    int spliced_read_len = strlen(seq_str); //FIXME TODO: use CIGAR instead
1295 >                    int left_splice_pos = atoi(splice_toks[0].c_str()) - left + 1;
1296                      if(left_splice_pos > spliced_read_len) left_splice_pos = spliced_read_len;
1297 <                    int8_t right_splice_pos = spliced_read_len - left_splice_pos;
1297 >                    int right_splice_pos = spliced_read_len - left_splice_pos;
1298  
1299                      if (right_splice_pos <= 0 || left_splice_pos <= 0)
1300                        return false;
# Line 1225 | Line 1313
1313                      //atoi(toks[right_window_edge_field].c_str());
1314                      int gap_len = atoi(splice_toks[1].c_str()) - atoi(splice_toks[0].c_str()) - 1;
1315  
1228                    string junction_strand = toks[num_extra_toks + strand_field];
1229                    if (!(junction_strand == "rev" || junction_strand == "fwd"))
1230                        // || !(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                      }
1237
1316                      int mismatches_in_anchor = 0;
1317 <                    for (int i=0;i<num_mismatches;i++) {
1318 <                                spl_num_mismatches++;
1319 <                                int mismatch_pos = mismatches[i];
1320 <                                if (((sam_flag & 0x0010) == 0 && abs(mismatch_pos - left_splice_pos) < (int)min_anchor_len) ||
1321 <                                  ((sam_flag & 0x0010) != 0 &&
1317 >                    for (size_t i=0;i<mismatches.size();++i) {
1318 >                                  spl_num_mismatches++;
1319 >                                  int mismatch_pos = mismatches[i];
1320 >                                  if (((sam_flag & 0x0010) == 0 && abs(mismatch_pos - left_splice_pos) < (int)min_anchor_len) ||
1321 >                                    ((sam_flag & 0x0010) != 0 &&
1322                                                    abs(((int)spliced_read_len - left_splice_pos + 1) - mismatch_pos)) < (int)min_anchor_len)
1323                                            mismatches_in_anchor++;
1324 <                                 }
1324 >                                  }
1325  
1326                      // FIXME: we probably should exclude these hits somewhere, but this
1327                      // isn't the right place
1328                      vector<CigarOp> splcigar;
1329                      CigarOpCode opcode=(toks[num_extra_toks + junction_type_field] == "del")? DEL : REF_SKIP ;
1330 <                    spliceCigar(splcigar, samcigar, left_splice_pos, gap_len, right_splice_pos, opcode);
1330 >                    //TODO: FIXME: use relative offsets for this call instead, like for the INS case
1331 >                    spl_num_mismatches=spliceCigar(splcigar, samcigar, mismatches, left, left_splice_pos, gap_len, opcode);
1332              /*
1333                      splcigar.push_back(CigarOp(MATCH, left_splice_pos));
1334                      if(toks[num_extra_toks + junction_type_field] == "del"){

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines