ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/tophat_cpp/bwt_map.cpp
(Generate patch)
# Line 1120 | Line 1120
1120   cigar.push_back(op);
1121   }
1122  
1123 < int spliceCigar(vector<CigarOp>& splcigar, vector<CigarOp>& cigar, vector<bool> mismatches,
1124 <                int &left, int spl_start, int spl_len, CigarOpCode spl_code) {
1125 < //merge the original 'cigar' with the new insert/gap operation
1126 < //at position spl_start and place the result into splcigar;
1127 < //TODO: ideally this should also get and rebuild the MD string (alignment mismatches)
1128 < int spl_mismatches=0;
1129 < //return value: mismatches in the insert region for INS case,
1130 < //or number of mismatches in the anchor region
1131 < //return -1 if somehow the hit seems bad
1132 <
1133 < //these offsets are relative to the beginning of alignment on reference
1134 < int spl_ofs=abs(spl_start-left); //relative position of splice op
1135 < int spl_ofs_end=spl_ofs; //relative position of first ref base AFTER splice op
1136 < CigarOp gapop(spl_code, spl_len); //for DEL, REF_SKIP, FUSIONS
1137 < if (spl_code==INS)
1138 <     spl_ofs_end += spl_len;
1139 < int ref_ofs=0; //working offset on reference
1140 < int read_ofs=0; //working offset on the read, relative to the leftmost aligned base
1141 < bool xfound=false;
1142 < //if (left<=spl_start+spl_len) {
1143 < if (spl_ofs_end>0) {
1144 <     int prev_opcode=0;
1145 <     int prev_oplen=0;
1146 <     for (size_t c = 0 ; c < cigar.size(); ++c)
1123 > bool spliceCigar(vector<CigarOp>& splcigar, vector<CigarOp>& cigar, vector<bool> mismatches,
1124 >                 int &left, int spl_start, int spl_len, CigarOpCode spl_code, int& spl_mismatches) {
1125 >  //merge the original 'cigar' with the new insert/gap operation
1126 >  //at position spl_start and place the result into splcigar;
1127 >  //TODO: ideally this should also get and rebuild the MD string (alignment mismatches)
1128 >
1129 >  //return value: mismatches in the insert region for INS case,
1130 >  //or number of mismatches in the anchor region
1131 >  //return -1 if somehow the hit seems bad
1132 >  
1133 >  //these offsets are relative to the beginning of alignment on reference
1134 >  int spl_ofs=spl_start-left; //relative position of splice op
1135 >  if (spl_code == FUSION_FF || spl_code == FUSION_FR || spl_code == FUSION_RF || spl_code == FUSION_RR)
1136 >    spl_ofs = abs(spl_ofs);
1137 >  int spl_ofs_end=spl_ofs; //relative position of first ref base AFTER splice op
1138 >  CigarOp gapop(spl_code, spl_len); //for DEL, REF_SKIP, FUSIONS
1139 >  if (spl_code==INS)
1140 >    spl_ofs_end += spl_len;
1141 >  
1142 >  int ref_ofs=0; //working offset on reference
1143 >  int read_ofs=0; //working offset on the read, relative to the leftmost aligned base
1144 >  bool xfound=false;
1145 >  //if (left<=spl_start+spl_len) {
1146 >  if (spl_ofs_end>0) {
1147 >    int prev_opcode=0;
1148 >    int prev_oplen=0;
1149 >    for (size_t c = 0 ; c < cigar.size(); ++c)
1150        {
1151          int prev_read_ofs=read_ofs;
1152          int cur_op_ofs=ref_ofs;
# Line 1184 | Line 1187
1187             }
1188  
1189          if (cur_op_ofs>=spl_ofs_end || ref_ofs<=spl_ofs) {
1190 <           if (cur_op_ofs==spl_ofs_end && spl_code!=INS)
1191 <              {
1192 <              xfound=true;
1193 <               //we have to insert the gap here first
1194 <              cigar_add(splcigar, gapop);
1195 <              //also, check
1196 <              }
1197 <           cigar_add(splcigar, cigar[c]);
1198 <           }
1190 >          if (cur_op_ofs==spl_ofs_end) {
1191 >            if (spl_code!=INS) {
1192 >              if (cur_opcode!=INS) {
1193 >                xfound=true;
1194 >                //we have to insert the gap here first
1195 >                cigar_add(splcigar, gapop);
1196 >                //also, check
1197 >              }
1198 >            }
1199 >          }
1200 >
1201 >          CigarOp op(cigar[c]);
1202 >          if (xfound)
1203 >            {
1204 >              if (spl_code == FUSION_FR || spl_code == FUSION_RR)
1205 >                {
1206 >                  if (op.opcode == MATCH)
1207 >                    op.opcode = mATCH;
1208 >                  else if (op.opcode == INS)
1209 >                    op.opcode = iNS;
1210 >                  else if (op.opcode == DEL)
1211 >                    op.opcode = dEL;
1212 >                  else if (op.opcode == REF_SKIP)
1213 >                    op.opcode = rEF_SKIP;
1214 >                }
1215 >            }
1216 >          else
1217 >            {
1218 >              if (spl_code == FUSION_RF || spl_code == FUSION_RR)
1219 >                {
1220 >                  if (op.opcode == MATCH)
1221 >                    op.opcode = mATCH;
1222 >                  else if (op.opcode == INS)
1223 >                    op.opcode = iNS;
1224 >                  else if (op.opcode == DEL)
1225 >                    op.opcode = dEL;
1226 >                  else if (op.opcode == REF_SKIP)
1227 >                    op.opcode = rEF_SKIP;
1228 >                }
1229 >            }
1230 >          
1231 >          cigar_add(splcigar, op);
1232 >        }
1233          else //if (ref_ofs>spl_ofs) {
1234             { //op intersection
1235             xfound=true;
# Line 1226 | Line 1263
1263                   CigarOp op(cigar[c]);
1264                   CigarOpCode opcode = op.opcode;
1265                   op.length=spl_ofs-cur_op_ofs;
1266 <                 if (gapop.opcode == FUSION_RF || gapop.opcode == FUSION_RR)
1266 >                 if (spl_code == FUSION_RF || spl_code == FUSION_RR)
1267                     {
1268                       if (opcode == MATCH)
1269                         op.opcode = mATCH;
# Line 1242 | Line 1279
1279                   cigar_add(splcigar, gapop);
1280  
1281                   op.opcode = opcode;
1282 <                 if (gapop.opcode == FUSION_FR || gapop.opcode == FUSION_RR)
1282 >                 if (spl_code == FUSION_FR || spl_code == FUSION_RR)
1283                     {
1284                       if (opcode == MATCH)
1285                         op.opcode = mATCH;
# Line 1277 | Line 1314
1314     //return spl_mismatches;
1315    // }
1316  
1317 < return spl_mismatches;
1317 >   return xfound;
1318   }
1319  
1320   bool SplicedSAMHitFactory::get_hit_from_buf(const char* orig_bwt_buf,
# Line 1420 | Line 1457
1457  
1458          vector<CigarOp> splcigar;
1459          //this also updates left to the adjusted genomic coordinates
1460 <        int spl_num_mismatches=spliceCigar(splcigar, samcigar, mismatches,
1461 <                                           left, left_splice_pos+1, insertedSequence.length(), INS);
1460 >        int spl_num_mismatches=0;
1461 >        bool overlapped = spliceCigar(splcigar, samcigar, mismatches,
1462 >                                      left, left_splice_pos+1, insertedSequence.length(),
1463 >                                      INS, spl_num_mismatches);
1464 >        
1465 >        if (!overlapped)
1466 >          return false;
1467  
1468          if (spl_num_mismatches<0) return false;
1469          num_mismatches-=spl_num_mismatches;
# Line 1551 | Line 1593
1593  
1594          CigarOpCode opcode=(toks[num_extra_toks + junction_type_field] == "del")? DEL : REF_SKIP;
1595          
1596 <        int spl_num_mismatches=spliceCigar(splcigar, samcigar, mismatches, left,
1597 <                                           left_splice_pos+1, gap_len, opcode);
1596 >        int spl_num_mismatches=0;
1597 >        bool overlapped = spliceCigar(splcigar, samcigar, mismatches, left,
1598 >                                      left_splice_pos+1, gap_len, opcode, spl_num_mismatches);
1599 >
1600 >        if (!overlapped)
1601 >          return false;
1602  
1603          if (spl_num_mismatches<0) // || spl_num_mismatches>max_anchor_mismatches)
1604            return false;
# Line 1593 | Line 1639
1639  
1640  
1641   bool BAMHitFactory::get_hit_from_buf(const char* orig_bwt_buf,
1642 <                                 BowtieHit& bh, bool strip_slash,
1643 <                                char* name_out, char* name_tags,
1644 <                                char* seq, char* qual) {
1642 >                                     BowtieHit& bh, bool strip_slash,
1643 >                                     char* name_out, char* name_tags,
1644 >                                     char* seq, char* qual) {
1645    if (_sam_header==NULL)
1646      err_die("Error: no SAM header when BAMHitFactory::get_hit_from_buf()!");
1647    const bam1_t* hit_buf = (const bam1_t*)orig_bwt_buf;
# Line 1813 | Line 1859
1859          int length = bam1_cigar(hit_buf)[i] >> BAM_CIGAR_SHIFT;
1860          if (length <= 0)
1861            {
1862 <            fprintf (stderr, "BAM error: CIGAR op has zero length\n");
1862 >            fprintf (stderr, "insert_id: %s - BAM error: CIGAR op has zero length\n", qname);
1863 >
1864 >            // daehwan - remove this
1865 >            exit(1);
1866 >            
1867              return false;
1868            }
1869          
# Line 2101 | Line 2151
2151              
2152              vector<CigarOp> splcigar;
2153              //this also updates left to the adjusted genomic coordinates
2154 <            int spl_num_mismatches=spliceCigar(splcigar, samcigar, mismatches,
2155 <                                               left, left_splice_pos+1, insertedSequence.length(), INS);
2154 >            int spl_num_mismatches=0;
2155 >            bool overlapped=spliceCigar(splcigar, samcigar, mismatches,
2156 >                                        left, left_splice_pos+1, insertedSequence.length(), INS, spl_num_mismatches);
2157 >            
2158 >            if (!overlapped)
2159 >              return false;
2160              
2161              if (spl_num_mismatches<0) return false;
2162              num_mismatches-=spl_num_mismatches;
# Line 2116 | Line 2170
2170                              num_mismatches,
2171                              0,
2172                              end);
2173 +
2174 +            // daehwan - remove this
2175 +            if (samcigar.size() > 1 && false)
2176 +              {
2177 +                cout << text_name << "\t" << left << endl
2178 +                     << print_cigar(samcigar) << endl
2179 +                     << print_cigar(splcigar) << endl;
2180 +              }
2181              
2182              return true;
2183            } //"ins"
# Line 2167 | Line 2229
2229                left_splice_pos -= 1;
2230              else
2231                left_splice_pos += 1;
2232 <            
2233 <            int spl_num_mismatches=spliceCigar(splcigar, samcigar, mismatches, left,
2234 <                                               left_splice_pos, gap_len, opcode);
2232 >
2233 >            int spl_num_mismatches=0;
2234 >            bool overlapped=spliceCigar(splcigar, samcigar, mismatches, left,
2235 >                                        left_splice_pos, gap_len, opcode, spl_num_mismatches);
2236 >
2237 >            if (!overlapped)
2238 >              return false;
2239              
2240              if (spl_num_mismatches<0) // || spl_num_mismatches>max_anchor_mismatches)
2241                return false;
# Line 2200 | Line 2266
2266                              spl_num_mismatches,
2267                              end);
2268  
2269 +            // daehwan - remove this
2270 +            if (samcigar.size() > 1 && false)
2271 +              {
2272 +                cout << text_name << "\t" << left << endl
2273 +                     << print_cigar(samcigar) << endl
2274 +                     << print_cigar(splcigar) << endl;
2275 +              }
2276 +
2277              return true;
2278            }
2279        } //parse splice data
# Line 2330 | Line 2404
2404          }
2405   }
2406  
2333 void print_hit(FILE* fout,
2334               const char* read_name,
2335               const BowtieHit& bh,
2336               const char* ref_name,
2337               const char* ref_name2,
2338               const char* sequence,
2339               const char* qualities)
2340 {
2341  string seq;
2342  string quals;
2343  if (sequence)
2344    {
2345      seq = sequence;
2346      quals = qualities;
2347      seq.resize(bh.read_len());
2348      quals.resize(bh.read_len());
2349    }
2350  else
2351    {
2352      seq = "*";
2353    }
2354  
2355  if (qualities)
2356    {
2357      quals = qualities;
2358      quals.resize(bh.read_len());
2359    }
2360  else
2361    {
2362      quals = "*";      
2363    }
2364  
2365  uint32_t sam_flag = 0;
2366  if (bh.antisense_align())
2367    {
2368      sam_flag |= 0x0010; // BAM_FREVERSE
2369    }
2370  
2371  int left = bh.left();
2372  
2373  uint32_t map_quality = 255;
2374  char cigar[256] = {0};
2375  string mate_ref_name = "*";
2376  uint32_t mate_pos = 0;
2377  uint32_t insert_size = 0;
2378  
2379  const vector<CigarOp>& bh_cigar = bh.cigar();
2380  
2381  /*
2382   * In addition to calculating the cigar string,
2383   * we need to figure out how many in/dels are in the
2384   * sequence, so that we can give the correct
2385   * value for the NM tag
2386   */
2387  int indel_distance = 0;
2388  CigarOpCode fusion_dir = FUSION_NOTHING;
2389  for (size_t c = 0; c < bh_cigar.size(); ++c)
2390    {
2391      const CigarOp& op = bh_cigar[c];
2392
2393      char ibuf[64];
2394      sprintf(ibuf, "%d", op.length);
2395      switch(op.opcode)
2396        {
2397        case MATCH:
2398        case mATCH:
2399          strcat(cigar, ibuf);
2400          if (bh_cigar[c].opcode == MATCH)
2401            strcat(cigar, "M");
2402          else
2403            strcat(cigar, "m");
2404          break;
2405        case INS:
2406        case iNS:
2407          strcat(cigar, ibuf);
2408          if (bh_cigar[c].opcode == INS)
2409            strcat(cigar, "I");
2410          else
2411            strcat(cigar, "i");
2412          indel_distance += bh_cigar[c].length;
2413          break;
2414        case DEL:
2415        case dEL:
2416          strcat(cigar, ibuf);
2417          if (bh_cigar[c].opcode == DEL)
2418            strcat(cigar, "D");
2419          else
2420            strcat(cigar, "d");
2421          indel_distance += bh_cigar[c].length;
2422          break;
2423        case REF_SKIP:
2424        case rEF_SKIP:
2425          strcat(cigar, ibuf);
2426          if (bh_cigar[c].opcode == REF_SKIP)
2427            strcat(cigar, "N");
2428          else
2429            strcat(cigar, "n");
2430          break;
2431        case FUSION_FF:
2432        case FUSION_FR:
2433        case FUSION_RF:
2434        case FUSION_RR:
2435          fusion_dir = op.opcode;
2436          sprintf(ibuf, "%d", bh_cigar[c].length + 1);
2437          strcat(cigar, ibuf);
2438          strcat(cigar, "F");
2439          break;
2440        default:
2441          break;
2442        }
2443    }
2444
2445  char cigar1[256] = {0}, cigar2[256] = {0};
2446  string left_seq, right_seq, left_qual, right_qual;
2447  int left1 = -1, left2 = -1;
2448  extract_partial_hits(bh, seq, quals,
2449                       cigar1, cigar2, left_seq, right_seq,
2450                       left_qual, right_qual, left1, left2);
2451
2452  
2453  bool containsSplice = false;
2454  for (vector<CigarOp>::const_iterator itr = bh.cigar().begin(); itr != bh.cigar().end(); ++itr)
2455    {
2456      if (itr->opcode == REF_SKIP || itr->opcode == rEF_SKIP)
2457        {
2458          containsSplice = true;
2459          break;
2460        }
2461    }
2462  
2463  if (fusion_dir != FUSION_NOTHING)
2464    {
2465      fprintf(fout,
2466              "%s\t%d\t%s\t%d\t%d\t%s\t%s\t%d\t%d\t%s\t%s\tNM:i:%d\t",
2467              read_name,
2468              sam_flag,
2469              ref_name,
2470              left1 + 1,
2471              map_quality,
2472              cigar1,
2473              mate_ref_name.c_str(),
2474              mate_pos,
2475              insert_size,
2476              left_seq.c_str(),
2477              left_qual.c_str(),
2478              bh.edit_dist() + indel_distance);
2479
2480      if (containsSplice)
2481        {
2482          fprintf(fout,
2483                  "XS:A:%c\t",
2484                  bh.antisense_splice() ? '-' : '+');
2485        }
2486
2487      fprintf(fout,
2488              "FR:Z:1 %s-%s %d %s %s %s\n",
2489              ref_name,
2490              ref_name2,
2491              left + 1,
2492              cigar,
2493              seq.c_str(),
2494              quals.c_str());
2495      
2496      fprintf(fout,
2497              "%s\t%d\t%s\t%d\t%d\t%s\t%s\t%d\t%d\t%s\t%s\tNM:i:%d\t",
2498              read_name,
2499              sam_flag,
2500              ref_name2,
2501              left2 + 1,
2502              map_quality,
2503              cigar2,
2504              mate_ref_name.c_str(),
2505              mate_pos,
2506              insert_size,
2507              right_seq.c_str(),
2508              right_qual.c_str(),
2509              bh.edit_dist() + indel_distance);
2510
2511      if (containsSplice)
2512        {
2513          fprintf(fout,
2514                  "XS:A:%c\t",
2515                  bh.antisense_splice() ? '-' : '+');
2516        }      
2517
2518      fprintf(fout,
2519              "FR:Z:2 %s-%s %d %s %s %s\n",
2520              ref_name,
2521              ref_name2,
2522              left + 1,
2523              cigar,
2524              seq.c_str(),
2525              quals.c_str());
2526    }
2527  else
2528    {
2529      fprintf(fout,
2530              "%s\t%d\t%s\t%d\t%d\t%s\t%s\t%d\t%d\t%s\t%s\tNM:i:%d\t",
2531              read_name,
2532              sam_flag,
2533              ref_name,
2534              left + 1,
2535              map_quality,
2536              cigar,
2537              mate_ref_name.c_str(),
2538              mate_pos,
2539              insert_size,
2540              seq.c_str(),
2541              quals.c_str(),
2542              bh.edit_dist() + indel_distance);
2543
2544      if (containsSplice)
2545        {
2546          fprintf(fout,
2547                  "XS:A:%c\t",
2548                  bh.antisense_splice() ? '-' : '+');
2549        }
2550
2551      fprintf(fout, "\n");
2552    }
2553 }
2554
2407   void print_bamhit(GBamWriter& wbam,
2408                    const char* read_name,
2409                    const BowtieHit& bh,
# Line 2559 | Line 2411
2411                    const char* ref_name2,
2412                    const char* sequence,
2413                    const char* qualities,
2414 <                  bool from_bowtie)
2414 >                  bool from_bowtie,
2415 >                  const vector<string>* extra_fields)
2416   {
2417    string seq;
2418    string quals;
# Line 2683 | Line 2536
2536      }
2537  
2538    vector<string> auxdata;
2539 +  if (extra_fields)
2540 +    auxdata.insert(auxdata.end(), extra_fields->begin(), extra_fields->end());
2541 +
2542    if (!sam_readgroup_id.empty())
2543      {
2544        string nm("RG:Z:");
# Line 2695 | Line 2551
2551    auxdata.push_back(nm);
2552    
2553    if (containsSplice) {
2554 <    nm="XS:A:";
2555 <    nm+=(char)(bh.antisense_splice() ? '-' : '+');
2556 <    auxdata.push_back(nm);
2554 >    // do not add more than once
2555 >    bool XS_found = false;
2556 >    for (size_t i = 0; i < auxdata.size(); ++i)
2557 >      {
2558 >        if (auxdata[i].substr(0, 2) == "XS")
2559 >          {
2560 >            XS_found = true;
2561 >            break;
2562 >          }
2563 >      }
2564 >
2565 >    if (!XS_found)
2566 >      {
2567 >        nm="XS:A:";
2568 >        nm+=(char)(bh.antisense_splice() ? '-' : '+');
2569 >        auxdata.push_back(nm);
2570 >      }
2571    }
2572    
2573    if (fusion_dir != FUSION_NOTHING)
# Line 2730 | Line 2600
2600                quals.c_str());
2601        auxdata.back() = XF;
2602      
2603 <      brec = wbam.new_record(read_name, sam_flag, ref_name, left2 + 1, map_quality,
2603 >      brec = wbam.new_record(read_name, sam_flag, ref_name2, left2 + 1, map_quality,
2604                               cigar2, mate_ref_name.c_str(), mate_pos,
2605                               insert_size, right_seq.c_str(), right_qual.c_str(), &auxdata);
2606        
# Line 3002 | Line 2872
2872   }
2873  
2874  
2875 < bool BowtieHit::check_editdist_consistency(const RefSequenceTable& rt)
2875 > bool BowtieHit::check_editdist_consistency(const RefSequenceTable& rt, bool bDebug)
2876   {
2877    RefSequenceTable::Sequence* ref_str1 = rt.get_seq(_ref_id);
2878    RefSequenceTable::Sequence* ref_str2 = rt.get_seq(_ref_id2);
2879    
2880    if (!ref_str1 || !ref_str2)
2881      return false;
2882 +
2883 +  if (bDebug)
2884 +    {
2885 +      cout << "check_editdist_consistency" << endl
2886 +           << "insert id: " << _insert_id << endl;
2887 +    }
2888    
2889    RefSequenceTable::Sequence* ref_str = ref_str1;
2890  
# Line 3026 | Line 2902
2902              seqan::Dna5String ref_seq = seqan::infix(*ref_str, pos_ref, pos_ref + cigar.length);
2903              for (size_t j = 0; j < cigar.length; ++j)
2904                {
2905 <                if (_seq[pos_seq] != ref_seq[j])
2905 >                seqan::Dna5 ref_nt = _seq[pos_seq];
2906 >                if (ref_nt != ref_seq[j])
2907                    ++mismatch;
2908 +
2909 +                if (bDebug)
2910 +                  cout << pos_seq << "\t" << ref_nt << " vs. " << ref_seq[j] << "\t" << mismatch << endl;
2911 +            
2912                  ++pos_seq;
2913                }
2914  
# Line 3042 | Line 2923
2923  
2924              for (size_t j = 0; j < cigar.length; ++j)
2925                {
2926 <                if (_seq[pos_seq] != ref_seq[j])
2926 >                seqan::Dna5 ref_nt = _seq[pos_seq];
2927 >                if (ref_nt != ref_seq[j])
2928                    ++mismatch;
2929 +
2930 +                if (bDebug)
2931 +                  cout << pos_seq << "\t" << ref_nt << " vs. " << ref_seq[j] << "\t" << mismatch << endl;
2932 +            
2933                  ++pos_seq;
2934                }
2935  
# Line 3092 | Line 2978
2978          }
2979      }
2980  
2981 +  if (bDebug)
2982 +    cout << "mismatch (real) vs. (calculated):" << mismatch << " vs. " << (int)_edit_dist << endl;
2983 +
2984    return mismatch == _edit_dist;
2985   }
2986 +
2987 + void bowtie_sam_extra(const BowtieHit& bh, const RefSequenceTable& rt, vector<string>& fields)
2988 + {
2989 +  RefSequenceTable::Sequence* ref_str1 = rt.get_seq(bh.ref_id());
2990 +  RefSequenceTable::Sequence* ref_str2 = rt.get_seq(bh.ref_id2());
2991 +  
2992 +  if (!ref_str1 || !ref_str2)
2993 +    return;
2994 +  
2995 +  RefSequenceTable::Sequence* ref_str = ref_str1;
2996 +
2997 +  size_t pos_seq = 0;
2998 +  size_t pos_mismatch = 0;
2999 +  size_t pos_ref = bh.left();
3000 +  size_t mismatch = 0;
3001 +  size_t num_gap_opens = 0;
3002 +  size_t num_gap_conts = 0;
3003 +  bool bSawFusion = false;
3004 +
3005 +  int AS_score = 0;
3006 +  
3007 +  const vector<CigarOp>& cigars = bh.cigar();
3008 +  const string& seq = bh.seq();
3009 +  const string& qual = bh.qual();
3010 +
3011 +  string AS = "AS:i:";
3012 +  string MD = "MD:Z:";
3013 +  
3014 +  for (size_t i = 0; i < cigars.size(); ++i)
3015 +    {
3016 +      CigarOp cigar = cigars[i];
3017 +      switch(cigar.opcode)
3018 +        {
3019 +        case MATCH:
3020 +        case mATCH:
3021 +          {
3022 +            seqan::Dna5String ref_seq;
3023 +            if (cigar.opcode == MATCH)
3024 +              {
3025 +                ref_seq = seqan::infix(*ref_str, pos_ref, pos_ref + cigar.length);
3026 +                pos_ref += cigar.length;
3027 +              }
3028 +            else
3029 +              {
3030 +                ref_seq = seqan::infix(*ref_str, pos_ref - cigar.length + 1, pos_ref + 1);
3031 +                seqan::convertInPlace(ref_seq, seqan::FunctorComplement<seqan::Dna>());
3032 +                seqan::reverseInPlace(ref_seq);
3033 +                pos_ref -= cigar.length;
3034 +              }
3035 +            
3036 +            for (size_t j = 0; j < cigar.length; ++j)
3037 +              {
3038 +                seqan::Dna5 ref_nt = ref_seq[j];
3039 +                if (seq[pos_seq] != ref_nt)
3040 +                  {
3041 +                    ++mismatch;
3042 +
3043 +                    if (pos_seq < qual.length())
3044 +                      {
3045 +                        float penalty = (bowtie2_max_penalty - bowtie2_min_penalty) * min((int)(qual[pos_seq] - '!'), 40) / 40.0;
3046 +                        AS_score -= (int)penalty;
3047 +                      }
3048 +
3049 +                    str_appendInt(MD, (int)pos_mismatch);
3050 +                    MD.push_back((char)ref_nt);
3051 +                    pos_mismatch = 0;
3052 +                  }
3053 +                else
3054 +                  ++pos_mismatch;
3055 +                
3056 +                ++pos_seq;
3057 +              }
3058 +          }
3059 +          break;
3060 +
3061 +        case INS:
3062 +        case iNS:
3063 +          {
3064 +            pos_seq += cigar.length;
3065 +
3066 +            AS_score -= bowtie2_read_gap_open;
3067 +            AS_score -= (int)(bowtie2_read_gap_cont * cigar.length);
3068 +
3069 +            num_gap_opens += 1;
3070 +            num_gap_conts += cigar.length;
3071 +          }
3072 +          break;
3073 +          
3074 +        case DEL:
3075 +        case dEL:
3076 +          {
3077 +            AS_score -= bowtie2_ref_gap_open;
3078 +            AS_score -= (int)(bowtie2_ref_gap_cont * cigar.length);
3079 +
3080 +            num_gap_opens += 1;
3081 +            num_gap_conts += cigar.length;
3082 +              
3083 +            seqan::Dna5String ref_seq;
3084 +            if (cigar.opcode == DEL)
3085 +              {
3086 +                ref_seq = seqan::infix(*ref_str, pos_ref, pos_ref + cigar.length);
3087 +                pos_ref += cigar.length;
3088 +              }
3089 +            else
3090 +              {
3091 +                ref_seq = seqan::infix(*ref_str, pos_ref - cigar.length + 1, pos_ref + 1);
3092 +                seqan::convertInPlace(ref_seq, seqan::FunctorComplement<seqan::Dna>());
3093 +                seqan::reverseInPlace(ref_seq);
3094 +                pos_ref -= cigar.length;
3095 +              }
3096 +
3097 +            str_appendInt(MD, (int)pos_mismatch);
3098 +            MD.push_back('^');
3099 +            for (size_t k = 0; k < length(ref_seq); ++k)
3100 +              MD.push_back((char)ref_seq[k]);
3101 +            
3102 +            pos_mismatch = 0;
3103 +          }
3104 +          break;
3105 +
3106 +        case REF_SKIP:
3107 +        case rEF_SKIP:
3108 +          {
3109 +            if (cigar.opcode == REF_SKIP)
3110 +              pos_ref += cigar.length;
3111 +            else
3112 +              pos_ref -= cigar.length;
3113 +          }
3114 +          break;
3115 +
3116 +        case FUSION_FF:
3117 +        case FUSION_FR:
3118 +        case FUSION_RF:
3119 +        case FUSION_RR:
3120 +          {
3121 +            // We don't allow a read spans more than two chromosomes.
3122 +            if (bSawFusion)
3123 +              return;
3124 +
3125 +            ref_str = ref_str2;  
3126 +            pos_ref = cigar.length;
3127 +
3128 +            bSawFusion = true;
3129 +          }
3130 +          break;
3131 +
3132 +        default:
3133 +          break;
3134 +        }
3135 +    }
3136 +
3137 +  str_appendInt(AS, AS_score);
3138 +  fields.push_back(AS);
3139 +
3140 +  string XM = "XM:i:";
3141 +  str_appendInt(XM, (int)mismatch);
3142 +  fields.push_back(XM);
3143 +
3144 +  string XO = "XO:i:";
3145 +  str_appendInt(XO, (int)num_gap_opens);
3146 +  fields.push_back(XO);
3147 +
3148 +  string XG = "XG:i:";
3149 +  str_appendInt(XG, (int)num_gap_conts);
3150 +  fields.push_back(XG);
3151 +
3152 +  str_appendInt(MD, (int)pos_mismatch);
3153 +  fields.push_back(MD);
3154 + }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines