ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/tophat_cpp/tophat_reports.cpp
(Generate patch)
# Line 88 | Line 88
88        // Is the new status better than the current best one?
89        if (best_grade < g)
90        {
91 <        best_hits.hits.clear();
91 >        if (!report_secondary_alignments)
92 >          best_hits.hits.clear();
93 >        
94          best_grade = g;
95          best_hits.hits.push_back(hits[i]);
96        }
# Line 158 | Line 160
160            // Is the new status better than the current best one?
161            if (best_grade < g)
162              {
163 <              left_best_hits.hits.clear();
164 <              right_best_hits.hits.clear();
163 >              if (!report_secondary_alignments)
164 >                {
165 >                  left_best_hits.hits.clear();
166 >                  right_best_hits.hits.clear();
167 >                }
168                
169                best_grade = g;
170                left_best_hits.hits.push_back(lh);
# Line 180 | Line 185
185   void add_auxData(vector<string>& auxdata, vector<string>& sam_toks,
186                   const RefSequenceTable& rt,const BowtieHit& bh, FragmentType insert_side,
187                   int num_hits, const BowtieHit* next_hit, int hitIndex) {
188 +  bool XS_found = false;
189 +  if (sam_toks.size()>11) {
190      
191 <    if (sam_toks.size()>11) {
192 <        for (size_t i=11;i<sam_toks.size();++i) {
193 <            if (sam_toks[i].find("NH:i:")==string::npos)
194 <                auxdata.push_back(sam_toks[i]);
195 <        }
196 <        }
197 <    string aux("NH:i:");
198 <    str_appendInt(aux, num_hits);
191 >    for (size_t i=11;i<sam_toks.size();++i) {
192 >      if (sam_toks[i].find("NH:i:")==string::npos &&
193 >          sam_toks[i].find("XS:i:")==string::npos)
194 >        auxdata.push_back(sam_toks[i]);
195 >
196 >      if (sam_toks[i].find("XS:A:")!=string::npos)
197 >        XS_found = true;
198 >    }
199 >  }
200 >  string aux("NH:i:");
201 >  str_appendInt(aux, num_hits);
202 >  auxdata.push_back(aux);
203 >  if (next_hit) {
204 >    const char* nh_ref_name = "=";
205 >    nh_ref_name = rt.get_name(next_hit->ref_id());
206 >    assert (nh_ref_name != NULL);
207 >    bool same_contig=(next_hit->ref_id()==bh.ref_id());
208 >    aux="CC:Z:";
209 >    aux+= (same_contig ? "=" : nh_ref_name);
210      auxdata.push_back(aux);
211 <    if (next_hit) {
212 <        const char* nh_ref_name = "=";
213 <        nh_ref_name = rt.get_name(next_hit->ref_id());
214 <        assert (nh_ref_name != NULL);
215 <        bool same_contig=(next_hit->ref_id()==bh.ref_id());
198 <        aux="CC:Z:";
199 <        aux+= (same_contig ? "=" : nh_ref_name);
200 <        auxdata.push_back(aux);
201 <        aux="CP:i:";
202 <        int nh_gpos=next_hit->left() + 1;
203 <        str_appendInt(aux, nh_gpos);
204 <        auxdata.push_back(aux);
205 <    } //has next_hit
211 >    aux="CP:i:";
212 >    int nh_gpos=next_hit->left() + 1;
213 >    str_appendInt(aux, nh_gpos);
214 >    auxdata.push_back(aux);
215 >  } //has next_hit
216      // FIXME: this code is still a bit brittle, because it contains no
217      // consistency check that the mates are on opposite strands - a current protocol
218      // requirement, and that the strand indicated by the alignment is consistent
219      // with the orientation of the splices (though that should be handled upstream).
220 +  if (!XS_found) {
221      const string xs_f("XS:A:+");
222      const string xs_r("XS:A:-");
223 <    if (bh.contiguous())  {
224 <        if (library_type == FR_FIRSTSTRAND) {
225 <            if (insert_side == FRAG_LEFT || insert_side == FRAG_UNPAIRED) {
226 <                if (bh.antisense_align())
227 <                    auxdata.push_back(xs_f);
228 <                else
229 <                    auxdata.push_back(xs_r);
230 <            }
231 <            else {
232 <                if (bh.antisense_align())
233 <                    auxdata.push_back(xs_r);
234 <                else
235 <                    auxdata.push_back(xs_f);
236 <            }
237 <        }
238 <        else if (library_type == FR_SECONDSTRAND)   {
239 <            if (insert_side == FRAG_LEFT || insert_side == FRAG_UNPAIRED){
240 <                if (bh.antisense_align())
241 <                    auxdata.push_back(xs_r);
242 <                else
243 <                    auxdata.push_back(xs_f);
244 <            }
245 <            else
246 <            {
247 <                if (bh.antisense_align())
248 <                    auxdata.push_back(xs_f);
249 <                else
250 <                    auxdata.push_back(xs_r);
251 <            }
252 <        }
253 <    } //bh.contiguous()
243 <    if (hitIndex >= 0)
223 >    if (library_type == FR_FIRSTSTRAND) {
224 >      if (insert_side == FRAG_LEFT || insert_side == FRAG_UNPAIRED) {
225 >        if (bh.antisense_align())
226 >          auxdata.push_back(xs_f);
227 >        else
228 >          auxdata.push_back(xs_r);
229 >      }
230 >      else {
231 >        if (bh.antisense_align())
232 >          auxdata.push_back(xs_r);
233 >        else
234 >          auxdata.push_back(xs_f);
235 >      }
236 >    }
237 >    else if (library_type == FR_SECONDSTRAND)   {
238 >      if (insert_side == FRAG_LEFT || insert_side == FRAG_UNPAIRED){
239 >        if (bh.antisense_align())
240 >          auxdata.push_back(xs_r);
241 >        else
242 >          auxdata.push_back(xs_f);
243 >      }
244 >      else
245 >        {
246 >          if (bh.antisense_align())
247 >            auxdata.push_back(xs_f);
248 >          else
249 >            auxdata.push_back(xs_r);
250 >        }
251 >    }
252 >  }
253 >  if (hitIndex >= 0)
254      {
255 <        string aux("HI:i:");
256 <        str_appendInt(aux, hitIndex);
257 <        auxdata.push_back(aux);
255 >      string aux("HI:i:");
256 >      str_appendInt(aux, hitIndex);
257 >      auxdata.push_back(aux);
258      }
259   }
260  
# Line 519 | Line 529
529  
530   void print_sam_for_single(const RefSequenceTable& rt,
531                            const HitsForRead& hits,
532 <                        const FragmentAlignmentGrade& grade,
533 <                        FragmentType frag_type,
534 <                        //FILE* reads_file,
535 <                        Read& read,
536 <                        GBamWriter& bam_writer,
537 <                        FILE* um_out//write unmapped reads here
528 <                        )
532 >                          const FragmentAlignmentGrade& grade,
533 >                          FragmentType frag_type,
534 >                          Read& read,
535 >                          GBamWriter& bam_writer,
536 >                          FILE* um_out //write unmapped reads here
537 >                          )
538   {
539      assert(!read.alt_name.empty());
540      lex_hit_sort s(rt, hits);
# Line 598 | Line 607
607        }
608      
609      bool got_left_read = left_reads_file.getRead(left_hits.insert_id, left_read,
610 <                                                 reads_format, false, left_um_out,
611 <                                                 false, begin_id, end_id);
610 >                                                 reads_format, false, begin_id, end_id,
611 >                                                 left_um_out, false);
612      
613      bool got_right_read = right_reads_file.getRead(right_hits.insert_id, right_read,
614 <                                                   reads_format, false, right_um_out,
615 <                                                   false, begin_id, end_id);
614 >                                                   reads_format, false, begin_id, end_id,
615 >                                                   right_um_out, false);
616      
617      if (left_hits.hits.size() == right_hits.hits.size())
618      {
# Line 1021 | Line 1030
1030              best_hits.insert_id = curr_left_obs_order;
1031              FragmentAlignmentGrade grade;
1032              bool got_read=left_reads_file.getRead(curr_left_obs_order, l_read,
1033 <                                                  reads_format, false, left_um_out.file,
1034 <                                                  false, begin_id, end_id);
1033 >                                                  reads_format, false, begin_id, end_id,
1034 >                                                  left_um_out.file, false);
1035              assert(got_read);
1036  
1037              if (right_reads_file.file()) {
1038                fprintf(left_um_out.file, "@%s #MAPPED#\n%s\n+\n%s\n", l_read.alt_name.c_str(),
1039                        l_read.seq.c_str(), l_read.qual.c_str());
1040                got_read=right_reads_file.getRead(curr_left_obs_order, r_read,
1041 <                                                reads_format, false, right_um_out.file,
1042 <                                                true, begin_id, end_id);
1041 >                                                reads_format, false, begin_id, end_id,
1042 >                                                right_um_out.file, true);
1043                assert(got_read);
1044              }
1045 <
1045 >    
1046              exclude_hits_on_filtered_junctions(*junctions, curr_left_hit_group);
1047  
1048              // Process hits for left singleton, select best alignments
# Line 1070 | Line 1079
1079              if (curr_right_obs_order >=  begin_id)
1080                {
1081                  bool got_read=right_reads_file.getRead(curr_right_obs_order, r_read,
1082 <                                                       reads_format, false, right_um_out.file,
1083 <                                                       false, begin_id, end_id);
1082 >                                                       reads_format, false, begin_id, end_id,
1083 >                                                       right_um_out.file, false);
1084  
1085                  assert(got_read);
1086  
1087                  fprintf(right_um_out.file, "@%s #MAPPED#\n%s\n+\n%s\n", r_read.alt_name.c_str(),
1088                          r_read.seq.c_str(), r_read.qual.c_str());
1089                  got_read=left_reads_file.getRead(curr_right_obs_order, l_read,
1090 <                                                 reads_format, false, left_um_out.file,
1091 <                                                 true, begin_id, end_id);
1090 >                                                 reads_format, false, begin_id, end_id,
1091 >                                                 left_um_out.file, true);
1092                  assert(got_read);
1093  
1094                  exclude_hits_on_filtered_junctions(*junctions, curr_right_hit_group);
# Line 1120 | Line 1129
1129                  //write it in the mapped file with the #MAPPED# flag
1130                  
1131                  bool got_read=right_reads_file.getRead(curr_left_obs_order, r_read,
1132 <                                                       reads_format, false, right_um_out.file,
1133 <                                                       false, begin_id, end_id);
1132 >                                                       reads_format, false, begin_id, end_id,
1133 >                                                       right_um_out.file, false);
1134                  assert(got_read);
1135                  fprintf(right_um_out.file, "@%s #MAPPED#\n%s\n+\n%s\n", r_read.alt_name.c_str(),
1136                          r_read.seq.c_str(), r_read.qual.c_str());
# Line 1152 | Line 1161
1161                  left_best_hits.insert_id = curr_left_obs_order;
1162                  //only left read mapped
1163                  bool got_read=left_reads_file.getRead(curr_left_obs_order, l_read,
1164 <                                                      reads_format, false, left_um_out.file,
1165 <                                                      false, begin_id, end_id);
1164 >                                                      reads_format, false, begin_id, end_id,
1165 >                                                      left_um_out.file, false);
1166                  assert(got_read);
1167                  fprintf(left_um_out.file, "@%s #MAPPED#\n%s\n+\n%s\n", l_read.alt_name.c_str(),
1168                          l_read.seq.c_str(), l_read.qual.c_str());
# Line 1224 | Line 1233
1233        } //while we still have unreported hits..
1234      //print the remaining unmapped reads at the end of each reads' stream
1235      left_reads_file.getRead(VMAXINT32, l_read,
1236 <                            reads_format, false, left_um_out.file,
1237 <                            false, begin_id, end_id);
1236 >                            reads_format, false, begin_id, end_id,
1237 >                            left_um_out.file, false);
1238      if (right_reads_file.file())
1239        right_reads_file.getRead(VMAXINT32, r_read,
1240 <                               reads_format, false, right_um_out.file,
1241 <                               false, begin_id, end_id);
1240 >                               reads_format, false, begin_id, end_id,
1241 >                               right_um_out.file, false);
1242  
1243      left_um_out.close();
1244      right_um_out.close();

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines