ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/tophat_cpp/fix_map_ordering.cpp
(Generate patch)
# Line 149 | Line 149
149                                          vector<pair<uint64_t, bam1_t*> >,
150                                          BamMapOrdering > map_pq;
151  
152 FILE* findex = NULL;
153 if (!index_outfile.empty()) {
154   findex = fopen(index_outfile.c_str(), "w");
155   if (findex == NULL)
156     err_die("Error: cannot create file %s\n", index_outfile.c_str());
157 }
158
152   tamFile fp=sam_open(fname.c_str());
153   bam1_t *b = bam_init1();
154   uint64_t last_id = 0;
162 uint64_t count = 0;
155   do {
156     bool new_record = (sam_read1(fp, header, b) >= 0);
157     if (new_record) {
# Line 175 | Line 167
167  
168       bool unmapped = (tb->core.flag & BAM_FUNMAP) != 0;
169       if (!unmapped) {
170 <       if (findex != NULL) {
171 <         if (count >= 1000 && t.first != last_id) {
180 <           int64_t offset = bam_writer.tell();
181 <           int block_offset = offset & 0xFFFF;
182 <          
183 <           // daehwan - this is a bug in bgzf.c in samtools
184 <           // I'll report this bug with a temporary solution, soon.
185 <           if (block_offset <= 0xF000) {
186 <             fprintf(findex, "%lu\t%ld\n", t.first, offset);
187 <             count = 0;
188 <           }
189 <         }
190 <
191 <         ++count;
192 <       }
193 <      
194 <       bam_writer.write(tb); //mapped
195 <     }
170 >       // mapped reads, write with index
171 >       bam_writer.write(tb, t.first);
172  
173 <     // In case of Bowtie2, some of the mapped reads against either transcriptome or genome
174 <     // may have low alignment scores due to gaps, in which case we will remap those.
175 <     // Later, we may have better alignments that usually involve splice junctions.
200 <     if (!unmapped) {
173 >       // In case of Bowtie2, some of the mapped reads against either transcriptome or genome
174 >       // may have low alignment scores due to gaps, in which case we will remap those.
175 >       // Later, we may have better alignments that usually involve splice junctions.
176         if (bowtie2 && t.first != last_id) {
177 <           uint8_t* ptr = bam_aux_get(tb, "AS");
178 <           if (ptr) {
179 <             int score = bam_aux2i(ptr);
180 <             if (score < bowtie2_min_score)
181 <               {
182 <                 unmapped = true;
183 <               }
184 <           }
177 >                   uint8_t* ptr = bam_aux_get(tb, "AS");
178 >                   if (ptr) {
179 >                         int score = bam_aux2i(ptr);
180 >                         if (score < bowtie2_min_score)
181 >                           {
182 >                                 unmapped = true;
183 >                           }
184 >                   }
185         }
211      
186         last_id = t.first;
187 <     }
187 >     } // mapped reads
188  
189 <     if (unmapped) { //unmapped
189 >     if (unmapped) { //unmapped reads or those who might be searched later
190         if (umbam!=NULL)
191 <         umbam->write(tb);
191 >         umbam->write(tb);
192       }
193       bam_destroy1(tb);
194       map_pq.pop();
# Line 224 | Line 198
198   bam_destroy1(b);
199   bam_header_destroy(header);
200  
227 if (findex != NULL)
228   fclose(findex);
201   }
202  
203   void driver_headerless(FILE* map_file, FILE* f_out)
# Line 244 | Line 216
216                  if (nl) *nl = 0;
217                  if (*bwt_buf == 0)
218                          continue;
219 <        /*
248 <                const char* bwt_fmt_str = "%s %c %s %d %s %s %d %s";
249 <                static const int buf_size = 2048;
250 <                char qname[buf_size];
251 <                char flagstr[buf_size];
252 <                int bwtf_ret = 0;
253 <                char refname[buf_size];
254 <                unsigned int text_offset;
255 <        char orientation;
256 <                char sequence[buf_size];
257 <                char qualities[buf_size];
258 <                unsigned int other_occs;
259 <                char mismatches[buf_size];
260 <                memset(mismatches, 0, sizeof(mismatches));
261 <                // Get a new record from the tab-delimited Bowtie map
262 <                bwtf_ret = sscanf(bwt_buf,
263 <                                                  bwt_fmt_str,
264 <                                                  qname,
265 <                                                  &orientation,
266 <                                                  refname,   // name of reference sequence
267 <                                                  &text_offset,
268 <                                                  sequence,
269 <                                                  qualities,
270 <                                                  &other_occs,
271 <                                                  mismatches);
272 <                
273 <                // If we didn't get enough fields, this record is bad, so skip it
274 <                if (bwtf_ret > 0 && bwtf_ret < 6)
275 <                {
276 <                        //fprintf(stderr, "Warning: found malformed record, skipping\n");
277 <                        continue;
278 <                }
279 <        */
280 <                TabSplitLine* l=new TabSplitLine(bwt_buf);
281 <                /* bwtf_ret = sscanf(bwt_buf, sam_fmt_str, qname, flagstr, refname);
282 <                if (qname[0]=='@' && qname[3]==0) { //header line, skip it
283 <                   continue;
284 <                   }
285 <                if (bwtf_ret>9 && bwtf_ret<3) {
286 <                   continue; //skip this line
287 <                   }
288 <                 */
219 >                TabSplitLine* l=new TabSplitLine(bwt_buf);
220                  if (l->tcount<10 || l->t[0][0]=='@') {
221                       delete l;
222                       continue;
# Line 363 | Line 294
294              }
295          else {
296                  //BAM output, we have SAM header
297 <                GBamWriter bam_writer(out_file_name.c_str(),sam_header.c_str());
297 >                GBamWriter bam_writer(out_file_name.c_str(),sam_header.c_str(), index_outfile);
298                  GBamWriter *unmapped_bam_writer=NULL;
299                  if (!out_unmapped_fname.empty())
300                      unmapped_bam_writer=new GBamWriter(out_unmapped_fname.c_str(),

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines