ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/tophat_cpp/bwt_map.cpp
(Generate patch)
# Line 773 | Line 773
773   }
774  
775   int getSAMmismatches(char* &buf, vector<CigarOp>& cigar,
776 <                  vector<int>& mismatches, int& sam_nm, bool& antisense_splice) {
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 809 | Line 809
809                    }
810                   while (isalpha(*p)) {
811                    p++;
812 <                  mismatches.push_back(bi);
812 >                  num_mismatches++;
813 >                  //mismatches.push_back(bi);
814 >                  mismatches[bi]=true;
815                    bi++;
816                    }
817                   if (*p=='^') { //reference deletion
# Line 820 | Line 822
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 850 | Line 852
852          break;
853                  }
854           }
855 <        return gspan;
855 >        return num_mismatches;
856   }
857  
858   bool SAMHitFactory::get_hit_from_buf(const char* orig_bwt_buf,
# Line 927 | Line 929
929          bool antisense_splice = false;
930          int sam_nm = 0; //the value of the NM tag (edit distance)
931          //int mismatches[1024];//array with mismatch positions on the read (0-based from the left aligned end of the read)
932 <        vector<int> mismatches;
933 <        getSAMmismatches(buf, cigar, mismatches, sam_nm, antisense_splice);
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 938 | Line 941
941                                  cigar,
942                                  sam_flag & 0x0010,
943                                  antisense_splice,
944 <                                mismatches.size(),
944 >                                num_mismatches,
945                                  0,
946                                  end);
947          }
# Line 951 | Line 954
954                                  cigar,
955                                  sam_flag & 0x0010,
956                                  false,
957 <                                mismatches.size(),
957 >                                num_mismatches,
958                                  0,
959                                  end);
960          }
# Line 966 | Line 969
969   cigar.push_back(op);
970   }
971  
972 < int spliceCigar(vector<CigarOp>& splcigar, vector<CigarOp>& cigar, vector<int> mismatches,
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 spl_start and place the result into splcigar;
976 < int num_mismatches=mismatches.size();
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
# Line 994 | Line 1002
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:
# Line 1006 | Line 1020
1020               break;
1021             //case HARD_CLIP:
1022             }
1023 +
1024          if (cur_op_ofs>=spl_ofs_end || ref_ofs<=spl_ofs) {
1025 <           if (spl_code!=INS && cur_op_ofs==spl_ofs_end)
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);
1029 >              cigar_add(splcigar, gapop);
1030 >              //also, check
1031                }
1032             cigar_add(splcigar, cigar[c]);
1033             }
# Line 1033 | Line 1050
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);
# Line 1045 | Line 1063
1063          prev_oplen=cur_oplen;
1064        } //for each cigar opcode
1065      } //intersection possible
1066 <  if (!xfound) {//no intersection found between splice event and alignment
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;
# Line 1055 | Line 1074
1074       //alignment ends before the splice event
1075       //nothing to do
1076      //  }
1077 <   }
1077 >   //return spl_mismatches;
1078 >  // }
1079  
1080 < //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;
1080 > return spl_mismatches;
1081   }
1082  
1083   bool SplicedSAMHitFactory::get_hit_from_buf(const char* orig_bwt_buf,
# Line 1135 | Line 1153
1153        return false;
1154          bool antisense_splice = false;
1155          int sam_nm = 0;
1156 <  vector<int> mismatches;
1157 <  getSAMmismatches(buf, samcigar, 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 1180 | Line 1200
1200        return false;
1201        }
1202  
1183    int spl_num_mismatches=0;
1184
1203                  //
1204                  // check for an insertion hit
1205                  //
# Line 1204 | Line 1222
1222  
1223        //TODO: implement this
1224        //this also updates left to the adjusted genomic coordinates
1225 <      spl_num_mismatches=spliceCigar(splcigar, samcigar, mismatches,
1226 <                        left, left_splice_pos, insertedSequence.length(), INS);
1225 >      int spl_num_mismatches=spliceCigar(splcigar, samcigar, mismatches,
1226 >                        left, left_splice_pos+1, insertedSequence.length(), INS);
1227 >      if (spl_num_mismatches<0) return false;
1228 >      num_mismatches-=spl_num_mismatches;
1229        /*
1230                          uint32_t right_splice_pos = left_splice_pos + 1;
1231  
# Line 1283 | Line 1303
1303                                          splcigar,
1304                                          sam_flag & 0x0010,
1305                                          junction_strand == "rev",
1306 <                                        spl_num_mismatches,
1306 >                                        num_mismatches,
1307                                          0,
1308                                          end);
1309                          return true;
1310                  } //"ins"
1311                  else //"del" or intron
1312                    {
1313 +                    // The 0-based position of the left edge of the alignment.
1314                      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                    int right_splice_pos = spliced_read_len - left_splice_pos;
1315  
1316 <                    if (right_splice_pos <= 0 || left_splice_pos <= 0)
1317 <                      return false;
1316 >                    // The 0-based position of the last genomic sequence before the deletion
1317 >              int left_splice_pos = atoi(splice_toks[0].c_str());
1318  
1319 +        int gap_len = atoi(splice_toks[1].c_str()) - left_splice_pos - 1;
1320 +        /*
1321                      if ((sam_flag & 0x0010) == 0) //######
1322                        {
1323 <            if (left_splice_pos + seg_offset < _anchor_length)
1323 >            if (left_splice_ofs + seg_offset < _anchor_length)
1324                return false;
1325            }
1326          else
1327            {
1328 <            if (right_splice_pos + seg_offset < _anchor_length)
1328 >            if (right_splice_ofs + seg_offset < _anchor_length)
1329                return false;
1330            }
1331 +          */
1332                      //uint32_t right = atoi(splice_toks[1].c_str()) + right_splice_pos;
1333                      //atoi(toks[right_window_edge_field].c_str());
1334 <                    int gap_len = atoi(splice_toks[1].c_str()) - atoi(splice_toks[0].c_str()) - 1;
1334 >
1335 >        /*
1336 >        //offset of deletion point, relative to the beginning of the alignment
1337 >        int left_splice_ofs = left_splice_pos-left+1;
1338  
1339                      int mismatches_in_anchor = 0;
1340                      for (size_t i=0;i<mismatches.size();++i) {
1341                                    spl_num_mismatches++;
1342                                    int mismatch_pos = mismatches[i];
1343 <                                  if (((sam_flag & 0x0010) == 0 && abs(mismatch_pos - left_splice_pos) < (int)min_anchor_len) ||
1343 >                                  if (((sam_flag & 0x0010) == 0 && abs(mismatch_pos - left_splice_ofs) < (int)min_anchor_len) ||
1344                                      ((sam_flag & 0x0010) != 0 &&
1345 <                                                  abs(((int)spliced_read_len - left_splice_pos + 1) - mismatch_pos)) < (int)min_anchor_len)
1345 >                                                  abs(((int)refspan - left_splice_ofs + 1) - mismatch_pos)) < (int)min_anchor_len)
1346                                            mismatches_in_anchor++;
1347                                    }
1348 <
1326 <                    // FIXME: we probably should exclude these hits somewhere, but this
1327 <                    // isn't the right place
1348 >         */
1349                      vector<CigarOp> splcigar;
1350                      CigarOpCode opcode=(toks[num_extra_toks + junction_type_field] == "del")? DEL : REF_SKIP ;
1351                      //TODO: FIXME: use relative offsets for this call instead, like for the INS case
1352 <                    spl_num_mismatches=spliceCigar(splcigar, samcigar, mismatches, left, left_splice_pos, gap_len, opcode);
1353 <            /*
1352 >                    int spl_num_mismatches=spliceCigar(splcigar, samcigar, mismatches, left,
1353 >                        left_splice_pos+1, gap_len, opcode);
1354 >                    if (spl_num_mismatches<0) // || spl_num_mismatches>max_anchor_mismatches)
1355 >                         return false;
1356 >        /*
1357                      splcigar.push_back(CigarOp(MATCH, left_splice_pos));
1358                      if(toks[num_extra_toks + junction_type_field] == "del"){
1359                        splcigar.push_back(CigarOp(DEL, gap_len));
# Line 1337 | Line 1361
1361                        splcigar.push_back(CigarOp(REF_SKIP, gap_len));
1362                      }
1363                      splcigar.push_back(CigarOp(MATCH, right_splice_pos));
1364 <            */
1364 >        */
1365                      bh = create_hit(name,
1366                                      contig,
1367                                      left,
1368                                      splcigar,
1369                                      (sam_flag & 0x0010),
1370                                      junction_strand == "rev",
1371 +                                    num_mismatches,
1372                                      spl_num_mismatches,
1348                                    mismatches_in_anchor,
1373                                      end);
1374                      return true;
1375                    }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines