1 |
/* |
2 |
* tophat_reports.cpp |
3 |
* TopHat |
4 |
* |
5 |
* Created by Cole Trapnell on 11/20/08. |
6 |
* Copyright 2008 Cole Trapnell. All rights reserved. |
7 |
* |
8 |
*/ |
9 |
|
10 |
#ifdef HAVE_CONFIG_H |
11 |
#include <config.h> |
12 |
#else |
13 |
#define PACKAGE_VERSION "INTERNAL" |
14 |
#define SVN_REVISION "XXX" |
15 |
#endif |
16 |
|
17 |
#include <cassert> |
18 |
#include <cstdio> |
19 |
#include <cstring> |
20 |
#include <vector> |
21 |
#include <string> |
22 |
#include <map> |
23 |
#include <algorithm> |
24 |
#include <set> |
25 |
#include <stdexcept> |
26 |
#include <iostream> |
27 |
#include <fstream> |
28 |
#include <seqan/sequence.h> |
29 |
#include <seqan/find.h> |
30 |
#include <seqan/file.h> |
31 |
#include <getopt.h> |
32 |
|
33 |
#include "common.h" |
34 |
#include "bwt_map.h" |
35 |
#include "junctions.h" |
36 |
#include "insertions.h" |
37 |
#include "deletions.h" |
38 |
#include "align_status.h" |
39 |
#include "fragments.h" |
40 |
#include "wiggles.h" |
41 |
#include "tokenize.h" |
42 |
#include "reads.h" |
43 |
|
44 |
#include "inserts.h" |
45 |
|
46 |
using namespace std; |
47 |
using namespace seqan; |
48 |
using std::set; |
49 |
|
50 |
void fragment_best_alignments(const HitsForRead& hits_for_read, |
51 |
FragmentAlignmentGrade& best_grade, |
52 |
HitsForRead& best_hits, |
53 |
const JunctionSet& gtf_junctions) |
54 |
{ |
55 |
const vector<BowtieHit>& hits = hits_for_read.hits; |
56 |
for (size_t i = 0; i < hits.size(); ++i) |
57 |
{ |
58 |
FragmentAlignmentGrade g(hits[i], gtf_junctions); |
59 |
// Is the new status better than the current best one? |
60 |
if (best_grade < g) |
61 |
{ |
62 |
best_hits.hits.clear(); |
63 |
best_grade = g; |
64 |
best_hits.hits.push_back(hits[i]); |
65 |
} |
66 |
else if (!(g < best_grade)) // is it just as good? |
67 |
{ |
68 |
best_grade.num_alignments++; |
69 |
best_hits.hits.push_back(hits[i]); |
70 |
} |
71 |
} |
72 |
} |
73 |
|
74 |
void insert_best_alignments(const HitsForRead& left_hits, |
75 |
const HitsForRead& right_hits, |
76 |
InsertAlignmentGrade& best_grade, |
77 |
HitsForRead& left_best_hits, |
78 |
HitsForRead& right_best_hits) |
79 |
{ |
80 |
// max mate inner distance (genomic) |
81 |
int min_mate_inner_dist = inner_dist_mean - inner_dist_std_dev; |
82 |
if (max_mate_inner_dist == -1) |
83 |
{ |
84 |
max_mate_inner_dist = inner_dist_mean + inner_dist_std_dev; |
85 |
} |
86 |
|
87 |
const vector<BowtieHit>& left = left_hits.hits; |
88 |
const vector<BowtieHit>& right = right_hits.hits; |
89 |
|
90 |
for (size_t i = 0; i < left.size(); ++i) |
91 |
{ |
92 |
const BowtieHit& lh = left[i]; |
93 |
for (size_t j = 0; j < right.size(); ++j) |
94 |
{ |
95 |
const BowtieHit& rh = right[j]; |
96 |
|
97 |
if (lh.ref_id() != rh.ref_id()) |
98 |
continue; |
99 |
|
100 |
InsertAlignmentGrade g(lh, rh, min_mate_inner_dist, max_mate_inner_dist); |
101 |
|
102 |
// Is the new status better than the current best one? |
103 |
if (best_grade < g) |
104 |
{ |
105 |
left_best_hits.hits.clear(); |
106 |
right_best_hits.hits.clear(); |
107 |
|
108 |
best_grade = g; |
109 |
left_best_hits.hits.push_back(lh); |
110 |
right_best_hits.hits.push_back(rh); |
111 |
} |
112 |
else if (!(g < best_grade)) |
113 |
{ |
114 |
left_best_hits.hits.push_back(lh); |
115 |
right_best_hits.hits.push_back(rh); |
116 |
best_grade.num_alignments++; |
117 |
} |
118 |
} |
119 |
} |
120 |
} |
121 |
|
122 |
enum FragmentType {FRAG_UNPAIRED, FRAG_LEFT, FRAG_RIGHT}; |
123 |
|
124 |
void add_auxData(vector<string>& auxdata, vector<string>& sam_toks, |
125 |
const RefSequenceTable& rt,const BowtieHit& bh, FragmentType insert_side, |
126 |
int num_hits, const BowtieHit* next_hit, int hitIndex) { |
127 |
|
128 |
if (sam_toks.size()>11) { |
129 |
for (size_t i=11;i<sam_toks.size();++i) { |
130 |
if (sam_toks[i].find("NH:i:")==string::npos) |
131 |
auxdata.push_back(sam_toks[i]); |
132 |
} |
133 |
} |
134 |
string aux("NH:i:"); |
135 |
str_appendInt(aux, num_hits); |
136 |
auxdata.push_back(aux); |
137 |
if (next_hit) { |
138 |
const char* nh_ref_name = "="; |
139 |
nh_ref_name = rt.get_name(next_hit->ref_id()); |
140 |
assert (nh_ref_name != NULL); |
141 |
bool same_contig=(next_hit->ref_id()==bh.ref_id()); |
142 |
aux="CC:Z:"; |
143 |
aux+= (same_contig ? "=" : nh_ref_name); |
144 |
auxdata.push_back(aux); |
145 |
aux="CP:i:"; |
146 |
int nh_gpos=next_hit->left() + 1; |
147 |
str_appendInt(aux, nh_gpos); |
148 |
auxdata.push_back(aux); |
149 |
} //has next_hit |
150 |
// FIXME: this code is still a bit brittle, because it contains no |
151 |
// consistency check that the mates are on opposite strands (a current protocol |
152 |
// requirement, and that the strand indicated by the alignment is consistent |
153 |
// with the orientation of the splices (though that should be handled upstream). |
154 |
const string xs_f("XS:A:+"); |
155 |
const string xs_r("XS:A:-"); |
156 |
if (bh.contiguous()) { |
157 |
if (library_type == FR_FIRSTSTRAND) { |
158 |
if (insert_side == FRAG_LEFT || insert_side == FRAG_UNPAIRED) { |
159 |
if (bh.antisense_align()) |
160 |
auxdata.push_back(xs_f); |
161 |
else |
162 |
auxdata.push_back(xs_r); |
163 |
} |
164 |
else { |
165 |
if (bh.antisense_align()) |
166 |
auxdata.push_back(xs_r); |
167 |
else |
168 |
auxdata.push_back(xs_f); |
169 |
} |
170 |
} |
171 |
else if (library_type == FR_SECONDSTRAND) { |
172 |
if (insert_side == FRAG_LEFT || insert_side == FRAG_UNPAIRED){ |
173 |
if (bh.antisense_align()) |
174 |
auxdata.push_back(xs_r); |
175 |
else |
176 |
auxdata.push_back(xs_f); |
177 |
} |
178 |
else |
179 |
{ |
180 |
if (bh.antisense_align()) |
181 |
auxdata.push_back(xs_f); |
182 |
else |
183 |
auxdata.push_back(xs_r); |
184 |
} |
185 |
} |
186 |
} //bh.contiguous() |
187 |
if (hitIndex >= 0) |
188 |
{ |
189 |
string aux("HI:i:"); |
190 |
str_appendInt(aux, hitIndex); |
191 |
auxdata.push_back(aux); |
192 |
} |
193 |
} |
194 |
|
195 |
bool rewrite_sam_record(GBamWriter& bam_writer, const RefSequenceTable& rt, |
196 |
const BowtieHit& bh, |
197 |
const char* bwt_buf, |
198 |
//char* rebuf, |
199 |
char* read_alt_name, |
200 |
const FragmentAlignmentGrade& grade, |
201 |
FragmentType insert_side, |
202 |
int num_hits, |
203 |
const BowtieHit* next_hit, |
204 |
bool primary, |
205 |
int hitIndex) |
206 |
{ |
207 |
// Rewrite this hit, filling in the alt name, mate mapping |
208 |
// and setting the pair flag |
209 |
vector<string> sam_toks; |
210 |
tokenize(bwt_buf, "\t", sam_toks); |
211 |
char* slash_pos = strrchr(read_alt_name,'/'); |
212 |
|
213 |
if (slash_pos) *slash_pos = 0; |
214 |
//*rebuf = 0; |
215 |
//read_alt_name as QNAME |
216 |
int flag=atoi(sam_toks[1].c_str()); //FLAG |
217 |
if (insert_side != FRAG_UNPAIRED) { |
218 |
//flag = atoi(sam_toks[1].c_str()); |
219 |
// mark this as a singleton mate |
220 |
flag |= 0x0001; |
221 |
if (insert_side == FRAG_LEFT) |
222 |
flag |= 0x0040; |
223 |
else if (insert_side == FRAG_RIGHT) |
224 |
flag |= 0x0080; |
225 |
flag |= 0x0008; |
226 |
//char flag_buf[64]; |
227 |
//sprintf(flag_buf, "%d", flag); |
228 |
//sam_toks[t] = flag_buf; |
229 |
} |
230 |
if (!primary) |
231 |
flag |= 0x100; |
232 |
|
233 |
int gpos=isdigit(sam_toks[3][0]) ? atoi(sam_toks[3].c_str()) : 0; |
234 |
int mapQ=255; |
235 |
if (grade.num_alignments > 1) { |
236 |
double err_prob = 1 - (1.0 / grade.num_alignments); |
237 |
mapQ = (int)(-10.0 * log(err_prob) / log(10.0)); |
238 |
} |
239 |
vector<string> auxdata; |
240 |
add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex); |
241 |
int tlen =atoi(sam_toks[8].c_str()); //TLEN |
242 |
int mate_pos=atoi(sam_toks[7].c_str()); |
243 |
GBamRecord* bamrec=bam_writer.new_record(read_alt_name, flag, sam_toks[2].c_str(), gpos, mapQ, |
244 |
sam_toks[5].c_str(), sam_toks[6].c_str(), mate_pos, |
245 |
tlen, sam_toks[9].c_str(), sam_toks[10].c_str(), &auxdata); |
246 |
bam_writer.write(bamrec); |
247 |
delete bamrec; |
248 |
return true; |
249 |
} |
250 |
|
251 |
bool rewrite_sam_record(GBamWriter& bam_writer, const RefSequenceTable& rt, |
252 |
const BowtieHit& bh, |
253 |
const char* bwt_buf, |
254 |
//char* rebuf, |
255 |
char* read_alt_name, |
256 |
const InsertAlignmentGrade& grade, |
257 |
FragmentType insert_side, |
258 |
const BowtieHit* partner, |
259 |
int num_hits, |
260 |
const BowtieHit* next_hit, |
261 |
bool primary, |
262 |
int hitIndex) |
263 |
{ |
264 |
// Rewrite this hit, filling in the alt name, mate mapping |
265 |
// and setting the pair flag |
266 |
vector<string> sam_toks; |
267 |
tokenize(bwt_buf, "\t", sam_toks); |
268 |
char* slash_pos = strrchr(read_alt_name,'/'); |
269 |
if (slash_pos) *slash_pos = 0; |
270 |
//read_alt_name as QNAME |
271 |
int flag = atoi(sam_toks[1].c_str()); |
272 |
// 0x0010 (strand of query) is assumed to be set correctly |
273 |
// to begin with |
274 |
flag |= 0x0001; //it must be paired |
275 |
if (insert_side == FRAG_LEFT) |
276 |
flag |= 0x0040; |
277 |
else if (insert_side == FRAG_RIGHT) |
278 |
flag |= 0x0080; |
279 |
if (!primary) |
280 |
flag |= 0x100; |
281 |
|
282 |
int gpos=isdigit(sam_toks[3][0]) ? atoi(sam_toks[3].c_str()) : 0; |
283 |
int mapQ=255; |
284 |
if (grade.num_alignments > 1) { |
285 |
double err_prob = 1 - (1.0 / grade.num_alignments); |
286 |
mapQ = (int)(-10.0 * log(err_prob) / log(10.0)); |
287 |
} |
288 |
int tlen=0; //TLEN |
289 |
int mate_pos=atoi(sam_toks[7].c_str()); |
290 |
if (partner) { |
291 |
if (partner->ref_id()==bh.ref_id()) { |
292 |
sam_toks[6] = "="; //same chromosome |
293 |
//TLEN: |
294 |
tlen = bh.left() < partner->left() ? partner->right() - bh.left() : |
295 |
partner->left() - bh.right(); |
296 |
} |
297 |
|
298 |
else { //partner on different chromosome/contig |
299 |
sam_toks[6] = rt.get_name(partner->ref_id()); |
300 |
if (sam_toks[6].empty()) { |
301 |
//FIXME -- this should never happen |
302 |
sam_toks[6] = "="; |
303 |
fprintf(stderr, "Warning: partner ref_id %d has no entry in ref table?\n", partner->ref_id()); |
304 |
} |
305 |
} |
306 |
mate_pos = partner->left() + 1; |
307 |
if (grade.happy()) |
308 |
flag |= 0x0002; |
309 |
if (partner->antisense_align()) |
310 |
flag |= 0x0020; |
311 |
} |
312 |
else { |
313 |
sam_toks[6] = "*"; |
314 |
mate_pos = 0; |
315 |
flag |= 0x0008; |
316 |
} |
317 |
vector<string> auxdata; |
318 |
add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex); |
319 |
GBamRecord* bamrec=bam_writer.new_record(read_alt_name, flag, sam_toks[2].c_str(), gpos, mapQ, |
320 |
sam_toks[5].c_str(), sam_toks[6].c_str(), mate_pos, |
321 |
tlen, sam_toks[9].c_str(), sam_toks[10].c_str(), &auxdata); |
322 |
bam_writer.write(bamrec); |
323 |
delete bamrec; |
324 |
return true; |
325 |
} |
326 |
|
327 |
struct lex_hit_sort |
328 |
{ |
329 |
lex_hit_sort(const RefSequenceTable& rt, const HitsForRead& hits) |
330 |
: _rt(rt), _hits(hits) |
331 |
{} |
332 |
|
333 |
bool operator()(const uint32_t& l, const uint32_t& r) const |
334 |
{ |
335 |
const BowtieHit& lhs = _hits.hits[l]; |
336 |
const BowtieHit& rhs = _hits.hits[r]; |
337 |
|
338 |
uint32_t l_id = lhs.ref_id(); |
339 |
uint32_t r_id = rhs.ref_id(); |
340 |
if (l_id != r_id) |
341 |
{ |
342 |
return (strcmp(_rt.get_name(lhs.ref_id()), _rt.get_name(rhs.ref_id())) < 0); |
343 |
} |
344 |
return lhs.left() < rhs.left(); |
345 |
} |
346 |
|
347 |
const RefSequenceTable& _rt; |
348 |
const HitsForRead& _hits; |
349 |
}; |
350 |
|
351 |
void print_sam_for_hits(const RefSequenceTable& rt, |
352 |
const HitsForRead& hits, |
353 |
const FragmentAlignmentGrade& grade, |
354 |
FragmentType frag_type, |
355 |
FILE* reads_file, |
356 |
GBamWriter& bam_writer) |
357 |
{ |
358 |
static const int buf_size = 2048; |
359 |
char read_name[buf_size]; |
360 |
char read_seq[buf_size]; |
361 |
char read_alt_name[buf_size]; |
362 |
char read_quals[buf_size]; |
363 |
|
364 |
//char rebuf[buf_size]; |
365 |
|
366 |
lex_hit_sort s(rt, hits); |
367 |
vector<uint32_t> index_vector; |
368 |
for (size_t i = 0; i < hits.hits.size(); ++i) |
369 |
index_vector.push_back(i); |
370 |
|
371 |
sort(index_vector.begin(), index_vector.end(), s); |
372 |
|
373 |
bool got_read = get_read_from_stream(hits.insert_id, |
374 |
reads_file, |
375 |
reads_format, |
376 |
false, |
377 |
read_name, |
378 |
read_seq, |
379 |
read_alt_name, |
380 |
read_quals); |
381 |
|
382 |
assert (got_read); |
383 |
|
384 |
size_t primaryHit = (hits.hits.empty()? 0: random() % hits.hits.size()); |
385 |
bool multipleHits = (hits.hits.size() > 1); |
386 |
for (size_t i = 0; i < hits.hits.size(); ++i) |
387 |
{ |
388 |
bool primary = (i == primaryHit); |
389 |
size_t index = index_vector[i]; |
390 |
const BowtieHit& bh = hits.hits[index]; |
391 |
rewrite_sam_record(bam_writer, rt, |
392 |
bh, |
393 |
bh.hitfile_rec().c_str(), |
394 |
read_alt_name, |
395 |
grade, |
396 |
frag_type, |
397 |
hits.hits.size(), |
398 |
(i < hits.hits.size()-1) ? &(hits.hits[index_vector[i+1]]) : NULL, |
399 |
primary, |
400 |
(multipleHits? i: -1)); |
401 |
} |
402 |
} |
403 |
|
404 |
void print_sam_for_hits(const RefSequenceTable& rt, |
405 |
const HitsForRead& left_hits, |
406 |
const HitsForRead& right_hits, |
407 |
const InsertAlignmentGrade& grade, |
408 |
FILE* left_reads_file, |
409 |
FILE* right_reads_file, |
410 |
GBamWriter& bam_writer) |
411 |
{ |
412 |
assert (left_hits.insert_id == right_hits.insert_id); |
413 |
|
414 |
static const int buf_size = 2048; |
415 |
char left_read_name[buf_size]; |
416 |
char left_read_seq[buf_size]; |
417 |
char left_read_alt_name[buf_size]; |
418 |
char left_read_quals[buf_size]; |
419 |
//char left_rebuf[buf_size]; |
420 |
|
421 |
char right_read_name[buf_size]; |
422 |
char right_read_seq[buf_size]; |
423 |
char right_read_alt_name[buf_size]; |
424 |
char right_read_quals[buf_size]; |
425 |
//char right_rebuf[buf_size]; |
426 |
|
427 |
assert (left_hits.hits.size() == right_hits.hits.size() || |
428 |
(left_hits.hits.empty() || right_hits.hits.empty())); |
429 |
|
430 |
size_t primaryHit = 0; |
431 |
vector<uint32_t> index_vector; |
432 |
if(right_hits.hits.size() > 0) |
433 |
{ |
434 |
lex_hit_sort s(rt, right_hits); |
435 |
for (size_t i = 0; i < right_hits.hits.size(); ++i) |
436 |
index_vector.push_back(i); |
437 |
|
438 |
sort(index_vector.begin(), index_vector.end(), s); |
439 |
primaryHit = random() % right_hits.hits.size(); |
440 |
} |
441 |
else if (left_hits.hits.size() > 0) |
442 |
{ |
443 |
lex_hit_sort s(rt, left_hits); |
444 |
for (size_t i = 0; i < left_hits.hits.size(); ++i) |
445 |
index_vector.push_back(i); |
446 |
|
447 |
sort(index_vector.begin(), index_vector.end(), s); |
448 |
primaryHit = random() % left_hits.hits.size(); |
449 |
} |
450 |
|
451 |
bool got_left_read = get_read_from_stream(left_hits.insert_id, |
452 |
left_reads_file, |
453 |
reads_format, |
454 |
false, |
455 |
left_read_name, |
456 |
left_read_seq, |
457 |
left_read_alt_name, |
458 |
left_read_quals); |
459 |
|
460 |
bool got_right_read = get_read_from_stream(right_hits.insert_id, |
461 |
right_reads_file, |
462 |
reads_format, |
463 |
false, |
464 |
right_read_name, |
465 |
right_read_seq, |
466 |
right_read_alt_name, |
467 |
right_read_quals); |
468 |
|
469 |
|
470 |
|
471 |
if (left_hits.hits.size() == right_hits.hits.size()) |
472 |
{ |
473 |
assert (got_left_read && got_right_read); |
474 |
bool multipleHits = (left_hits.hits.size() > 1); |
475 |
for (size_t i = 0; i < right_hits.hits.size(); ++i) |
476 |
{ |
477 |
bool primary = (i == primaryHit); |
478 |
size_t index = index_vector[i]; |
479 |
const BowtieHit& right_bh = right_hits.hits[index]; |
480 |
const BowtieHit& left_bh = left_hits.hits[index]; |
481 |
|
482 |
rewrite_sam_record(bam_writer, rt, |
483 |
right_bh, |
484 |
right_bh.hitfile_rec().c_str(), |
485 |
right_read_alt_name, |
486 |
grade, |
487 |
FRAG_RIGHT, |
488 |
&left_bh, |
489 |
right_hits.hits.size(), |
490 |
(i < right_hits.hits.size() - 1) ? &(right_hits.hits[index_vector[i+1]]) : NULL, |
491 |
primary, |
492 |
(multipleHits? i: -1)); |
493 |
rewrite_sam_record(bam_writer, rt, |
494 |
left_bh, |
495 |
left_bh.hitfile_rec().c_str(), |
496 |
left_read_alt_name, |
497 |
grade, |
498 |
FRAG_LEFT, |
499 |
&right_bh, |
500 |
left_hits.hits.size(), |
501 |
(i < left_hits.hits.size() - 1) ? &(left_hits.hits[index_vector[i+1]]) : NULL, |
502 |
primary, |
503 |
(multipleHits? i: -1)); |
504 |
} |
505 |
} |
506 |
else if (left_hits.hits.empty()) |
507 |
{ |
508 |
for (size_t i = 0; i < right_hits.hits.size(); ++i) |
509 |
{ |
510 |
bool primary = (i == primaryHit); |
511 |
size_t index = index_vector[i]; |
512 |
const BowtieHit& bh = right_hits.hits[index]; |
513 |
|
514 |
rewrite_sam_record(bam_writer, rt, |
515 |
bh, |
516 |
bh.hitfile_rec().c_str(), |
517 |
right_read_alt_name, |
518 |
grade, |
519 |
FRAG_RIGHT, |
520 |
NULL, |
521 |
right_hits.hits.size(), |
522 |
(i < right_hits.hits.size() - 1) ? &(right_hits.hits[index_vector[i+1]]) : NULL, |
523 |
primary, |
524 |
-1); |
525 |
|
526 |
} |
527 |
} |
528 |
else if (right_hits.hits.empty()) |
529 |
{ |
530 |
for (size_t i = 0; i < left_hits.hits.size(); ++i) |
531 |
{ |
532 |
bool primary = (i == primaryHit); |
533 |
size_t index = index_vector[i]; |
534 |
const BowtieHit& bh = left_hits.hits[index]; |
535 |
rewrite_sam_record(bam_writer, rt, |
536 |
bh, |
537 |
bh.hitfile_rec().c_str(), |
538 |
left_read_alt_name, |
539 |
grade, |
540 |
FRAG_LEFT, |
541 |
NULL, |
542 |
left_hits.hits.size(), |
543 |
(i < left_hits.hits.size() - 1) ? &(left_hits.hits[index_vector[i+1]]) : NULL, |
544 |
primary, |
545 |
-1); |
546 |
|
547 |
} |
548 |
} |
549 |
else |
550 |
{ |
551 |
assert (false); |
552 |
} |
553 |
} |
554 |
|
555 |
/** |
556 |
* Given all of the hits for a particular read, update the read counts for insertions and deletions. |
557 |
* @param hits hits The alignments for a particular read |
558 |
* @param insertions Maps from an insertion to the number of supporting reads for that insertion |
559 |
* @param deletions Maps from a deletion to the number of supporting reads for that deletion |
560 |
*/ |
561 |
void update_insertions_and_deletions(const HitsForRead& hits, |
562 |
InsertionSet& insertions, |
563 |
DeletionSet& deletions) |
564 |
{ |
565 |
for (size_t i = 0; i < hits.hits.size(); ++i) |
566 |
{ |
567 |
const BowtieHit& bh = hits.hits[i]; |
568 |
insertions_from_alignment(bh, insertions); |
569 |
deletions_from_alignment(bh, deletions); |
570 |
} |
571 |
} |
572 |
|
573 |
|
574 |
|
575 |
void update_junctions(const HitsForRead& hits, |
576 |
JunctionSet& junctions) |
577 |
{ |
578 |
for (size_t i = 0; i < hits.hits.size(); ++i) |
579 |
{ |
580 |
const BowtieHit& bh = hits.hits[i]; |
581 |
junctions_from_alignment(bh, junctions); |
582 |
} |
583 |
} |
584 |
|
585 |
// Extracts junctions from all the SAM hits (based on REF_SKIPs) in the hit file |
586 |
// resets the stream when finished. |
587 |
void exclude_hits_on_filtered_junctions(const JunctionSet& junctions, |
588 |
HitsForRead& hits) |
589 |
{ |
590 |
HitsForRead remaining; |
591 |
remaining.insert_id = hits.insert_id; |
592 |
|
593 |
for (size_t i = 0; i < hits.hits.size(); ++i) |
594 |
{ |
595 |
BowtieHit& bh = hits.hits[i]; |
596 |
bool filter_hit = false; |
597 |
if (!bh.contiguous()) |
598 |
{ |
599 |
JunctionSet bh_juncs; |
600 |
junctions_from_alignment(bh, bh_juncs); |
601 |
for (JunctionSet::iterator itr = bh_juncs.begin(); |
602 |
itr != bh_juncs.end(); |
603 |
itr++) |
604 |
{ |
605 |
const Junction& j = itr->first; |
606 |
JunctionSet::const_iterator target = junctions.find(j); |
607 |
if (target == junctions.end() || !target->second.accepted) |
608 |
{ |
609 |
filter_hit = true; |
610 |
break; |
611 |
} |
612 |
} |
613 |
} |
614 |
if (!filter_hit) |
615 |
{ |
616 |
remaining.hits.push_back(bh); |
617 |
} |
618 |
else |
619 |
{ |
620 |
//fprintf(stderr, "Filtering hit\n"); |
621 |
} |
622 |
} |
623 |
hits = remaining; |
624 |
} |
625 |
|
626 |
void get_junctions_from_best_hits(HitStream& left_hs, |
627 |
HitStream& right_hs, |
628 |
ReadTable& it, |
629 |
JunctionSet& junctions, |
630 |
const JunctionSet& gtf_junctions) |
631 |
{ |
632 |
HitsForRead curr_left_hit_group; |
633 |
HitsForRead curr_right_hit_group; |
634 |
|
635 |
left_hs.next_read_hits(curr_left_hit_group); |
636 |
right_hs.next_read_hits(curr_right_hit_group); |
637 |
|
638 |
uint32_t curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id); |
639 |
uint32_t curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id); |
640 |
|
641 |
// While we still have unreported hits... |
642 |
while(curr_left_obs_order != 0xFFFFFFFF || |
643 |
curr_right_obs_order != 0xFFFFFFFF) |
644 |
{ |
645 |
// Chew up left singletons |
646 |
while (curr_left_obs_order < curr_right_obs_order && |
647 |
curr_left_obs_order != 0xFFFFFFFF) |
648 |
{ |
649 |
HitsForRead best_hits; |
650 |
best_hits.insert_id = curr_left_obs_order; |
651 |
FragmentAlignmentGrade grade; |
652 |
|
653 |
// Process hits for left singleton, select best alignments |
654 |
fragment_best_alignments(curr_left_hit_group, grade, best_hits, gtf_junctions); |
655 |
update_junctions(best_hits, junctions); |
656 |
|
657 |
// Get next hit group |
658 |
left_hs.next_read_hits(curr_left_hit_group); |
659 |
curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id); |
660 |
} |
661 |
|
662 |
// Chew up right singletons |
663 |
while (curr_left_obs_order > curr_right_obs_order && |
664 |
curr_right_obs_order != 0xFFFFFFFF) |
665 |
{ |
666 |
HitsForRead best_hits; |
667 |
best_hits.insert_id = curr_right_obs_order; |
668 |
FragmentAlignmentGrade grade; |
669 |
|
670 |
// Process hit for right singleton, select best alignments |
671 |
fragment_best_alignments(curr_right_hit_group,grade, best_hits, gtf_junctions); |
672 |
update_junctions(best_hits, junctions); |
673 |
|
674 |
// Get next hit group |
675 |
right_hs.next_read_hits(curr_right_hit_group); |
676 |
curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id); |
677 |
} |
678 |
|
679 |
// Since we have both left hits and right hits for this insert, |
680 |
// Find the best pairing and print both |
681 |
while (curr_left_obs_order == curr_right_obs_order && |
682 |
curr_left_obs_order != 0xFFFFFFFF && curr_right_obs_order != 0xFFFFFFFF) |
683 |
{ |
684 |
if (curr_left_hit_group.hits.empty()) |
685 |
{ |
686 |
HitsForRead right_best_hits; |
687 |
right_best_hits.insert_id = curr_right_obs_order; |
688 |
|
689 |
FragmentAlignmentGrade grade; |
690 |
fragment_best_alignments(curr_right_hit_group, grade, right_best_hits, gtf_junctions); |
691 |
update_junctions(right_best_hits, junctions); |
692 |
} |
693 |
else if (curr_right_hit_group.hits.empty()) |
694 |
{ |
695 |
HitsForRead left_best_hits; |
696 |
left_best_hits.insert_id = curr_left_obs_order; |
697 |
|
698 |
FragmentAlignmentGrade grade; |
699 |
// Process hits for left singleton, select best alignments |
700 |
fragment_best_alignments(curr_left_hit_group, grade, left_best_hits, gtf_junctions); |
701 |
update_junctions(left_best_hits, junctions); |
702 |
} |
703 |
else |
704 |
{ |
705 |
HitsForRead left_best_hits; |
706 |
HitsForRead right_best_hits; |
707 |
left_best_hits.insert_id = curr_left_obs_order; |
708 |
right_best_hits.insert_id = curr_right_obs_order; |
709 |
|
710 |
// daehwan - apply gtf_junctions here, too! |
711 |
|
712 |
InsertAlignmentGrade grade; |
713 |
insert_best_alignments(curr_left_hit_group, |
714 |
curr_right_hit_group, |
715 |
grade, |
716 |
left_best_hits, |
717 |
right_best_hits); |
718 |
|
719 |
update_junctions(left_best_hits, junctions); |
720 |
update_junctions(right_best_hits, junctions); |
721 |
} |
722 |
|
723 |
left_hs.next_read_hits(curr_left_hit_group); |
724 |
curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id); |
725 |
|
726 |
right_hs.next_read_hits(curr_right_hit_group); |
727 |
curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id); |
728 |
} |
729 |
} |
730 |
left_hs.reset(); |
731 |
right_hs.reset(); |
732 |
} |
733 |
|
734 |
|
735 |
void driver(GBamWriter& bam_writer, |
736 |
string& left_map_fname, |
737 |
FILE* left_reads, |
738 |
string& right_map_fname, |
739 |
FILE* right_reads, |
740 |
FILE* junctions_out, |
741 |
FILE* insertions_out, |
742 |
FILE* deletions_out) |
743 |
{ |
744 |
ReadTable it; |
745 |
RefSequenceTable rt(sam_header, true); |
746 |
srandom(1); |
747 |
JunctionSet gtf_junctions; |
748 |
if (gtf_juncs != "") |
749 |
{ |
750 |
char splice_buf[2048]; |
751 |
FILE* splice_coords = fopen(gtf_juncs.c_str(), "r"); |
752 |
if (splice_coords) |
753 |
{ |
754 |
while (fgets(splice_buf, sizeof(splice_buf), splice_coords)) |
755 |
{ |
756 |
char* nl = strrchr(splice_buf, '\n'); |
757 |
char* buf = splice_buf; |
758 |
if (nl) *nl = 0; |
759 |
|
760 |
char* ref_name = get_token((char**)&buf, "\t"); |
761 |
char* scan_left_coord = get_token((char**)&buf, "\t"); |
762 |
char* scan_right_coord = get_token((char**)&buf, "\t"); |
763 |
char* orientation = get_token((char**)&buf, "\t"); |
764 |
|
765 |
if (!scan_left_coord || !scan_right_coord || !orientation) |
766 |
{ |
767 |
fprintf(stderr,"Error: malformed splice coordinate record\n"); |
768 |
exit(1); |
769 |
} |
770 |
|
771 |
uint32_t ref_id = rt.get_id(ref_name, NULL, 0); |
772 |
uint32_t left_coord = atoi(scan_left_coord); |
773 |
uint32_t right_coord = atoi(scan_right_coord); |
774 |
bool antisense = *orientation == '-'; |
775 |
|
776 |
// add 1 to left_coord to meet BED format |
777 |
gtf_junctions.insert(make_pair<Junction, JunctionStats>(Junction(ref_id, left_coord + 1, right_coord, antisense), JunctionStats())); |
778 |
} |
779 |
} |
780 |
} |
781 |
|
782 |
BAMHitFactory hit_factory(it,rt); |
783 |
JunctionSet junctions; |
784 |
{ |
785 |
HitStream l_hs(left_map_fname, &hit_factory, false, true, true, true); |
786 |
HitStream r_hs(right_map_fname, &hit_factory, false, true, true, true); |
787 |
get_junctions_from_best_hits(l_hs, r_hs, it, junctions, gtf_junctions); |
788 |
//this resets the streams |
789 |
} |
790 |
|
791 |
HitStream left_hs(left_map_fname, &hit_factory, false, true, true, true); |
792 |
HitStream right_hs(right_map_fname, &hit_factory, false, true, true, true); |
793 |
|
794 |
size_t num_unfiltered_juncs = junctions.size(); |
795 |
fprintf(stderr, "Loaded %lu junctions\n", (long unsigned int) num_unfiltered_juncs); |
796 |
|
797 |
HitsForRead curr_left_hit_group; |
798 |
HitsForRead curr_right_hit_group; |
799 |
|
800 |
left_hs.next_read_hits(curr_left_hit_group); |
801 |
right_hs.next_read_hits(curr_right_hit_group); |
802 |
|
803 |
uint32_t curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id); |
804 |
uint32_t curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id); |
805 |
|
806 |
// Read hits, extract junctions, and toss the ones that arent strongly enough supported. |
807 |
filter_junctions(junctions, gtf_junctions); |
808 |
//size_t num_juncs_after_filter = junctions.size(); |
809 |
//fprintf(stderr, "Filtered %lu junctions\n", num_unfiltered_juncs - num_juncs_after_filter); |
810 |
|
811 |
size_t small_overhangs = 0; |
812 |
for (JunctionSet::iterator i = junctions.begin(); i != junctions.end(); ++i) |
813 |
{ |
814 |
if (i->second.accepted && |
815 |
(i->second.left_extent < min_anchor_len || i->second.right_extent < min_anchor_len)) |
816 |
{ |
817 |
small_overhangs++; |
818 |
} |
819 |
} |
820 |
|
821 |
if (small_overhangs >0) |
822 |
fprintf(stderr, "Warning: %lu small overhang junctions!\n", (long unsigned int)small_overhangs); |
823 |
|
824 |
JunctionSet final_junctions; // the junctions formed from best hits |
825 |
InsertionSet final_insertions; |
826 |
DeletionSet final_deletions; |
827 |
|
828 |
fprintf (stderr, "Reporting final accepted alignments..."); |
829 |
// While we still have unreported hits... |
830 |
while(curr_left_obs_order != 0xFFFFFFFF || |
831 |
curr_right_obs_order != 0xFFFFFFFF) |
832 |
{ |
833 |
// Chew up left singletons |
834 |
while (curr_left_obs_order < curr_right_obs_order && |
835 |
curr_left_obs_order != 0xFFFFFFFF) |
836 |
{ |
837 |
HitsForRead best_hits; |
838 |
best_hits.insert_id = curr_left_obs_order; |
839 |
FragmentAlignmentGrade grade; |
840 |
|
841 |
exclude_hits_on_filtered_junctions(junctions, curr_left_hit_group); |
842 |
|
843 |
// Process hits for left singleton, select best alignments |
844 |
fragment_best_alignments(curr_left_hit_group, grade, best_hits, gtf_junctions); |
845 |
|
846 |
if (best_hits.hits.size() <= max_multihits) |
847 |
{ |
848 |
update_junctions(best_hits, final_junctions); |
849 |
update_insertions_and_deletions(best_hits, final_insertions, final_deletions); |
850 |
|
851 |
print_sam_for_hits(rt, |
852 |
best_hits, |
853 |
grade, |
854 |
(right_map_fname.empty() ? FRAG_UNPAIRED : FRAG_LEFT), |
855 |
left_reads, |
856 |
bam_writer); |
857 |
} |
858 |
|
859 |
// Get next hit group |
860 |
left_hs.next_read_hits(curr_left_hit_group); |
861 |
curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id); |
862 |
} //left singletons |
863 |
|
864 |
// Chew up right singletons |
865 |
while (curr_left_obs_order > curr_right_obs_order && |
866 |
curr_right_obs_order != 0xFFFFFFFF) |
867 |
{ |
868 |
HitsForRead best_hits; |
869 |
best_hits.insert_id = curr_right_obs_order; |
870 |
FragmentAlignmentGrade grade; |
871 |
|
872 |
exclude_hits_on_filtered_junctions(junctions, curr_right_hit_group); |
873 |
|
874 |
// Process hit for right singleton, select best alignments |
875 |
fragment_best_alignments(curr_right_hit_group, grade, best_hits, gtf_junctions); |
876 |
|
877 |
if (best_hits.hits.size() <= max_multihits) |
878 |
{ |
879 |
update_junctions(best_hits, final_junctions); |
880 |
update_insertions_and_deletions(best_hits, final_insertions, final_deletions); |
881 |
|
882 |
print_sam_for_hits(rt, |
883 |
best_hits, |
884 |
grade, |
885 |
FRAG_RIGHT, |
886 |
right_reads, |
887 |
bam_writer); |
888 |
} |
889 |
|
890 |
// Get next hit group |
891 |
right_hs.next_read_hits(curr_right_hit_group); |
892 |
curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id); |
893 |
} |
894 |
|
895 |
// Since we have both left hits and right hits for this insert, |
896 |
// Find the best pairing and print both |
897 |
while (curr_left_obs_order == curr_right_obs_order && |
898 |
curr_left_obs_order != 0xFFFFFFFF && curr_right_obs_order != 0xFFFFFFFF) |
899 |
{ |
900 |
exclude_hits_on_filtered_junctions(junctions, curr_left_hit_group); |
901 |
exclude_hits_on_filtered_junctions(junctions, curr_right_hit_group); |
902 |
|
903 |
if (curr_left_hit_group.hits.empty()) |
904 |
{ |
905 |
HitsForRead right_best_hits; |
906 |
right_best_hits.insert_id = curr_right_obs_order; |
907 |
|
908 |
FragmentAlignmentGrade grade; |
909 |
fragment_best_alignments(curr_right_hit_group, grade, right_best_hits, gtf_junctions); |
910 |
|
911 |
if (right_best_hits.hits.size() <= max_multihits) |
912 |
{ |
913 |
update_junctions(right_best_hits, final_junctions); |
914 |
update_insertions_and_deletions(right_best_hits, final_insertions, final_deletions); |
915 |
|
916 |
print_sam_for_hits(rt, |
917 |
right_best_hits, |
918 |
grade, |
919 |
FRAG_RIGHT, |
920 |
right_reads, |
921 |
bam_writer); |
922 |
} |
923 |
} |
924 |
else if (curr_right_hit_group.hits.empty()) |
925 |
{ |
926 |
HitsForRead left_best_hits; |
927 |
left_best_hits.insert_id = curr_left_obs_order; |
928 |
|
929 |
FragmentAlignmentGrade grade; |
930 |
// Process hits for left singleton, select best alignments |
931 |
fragment_best_alignments(curr_left_hit_group, grade, left_best_hits, gtf_junctions); |
932 |
|
933 |
if (left_best_hits.hits.size() <= max_multihits) |
934 |
{ |
935 |
update_junctions(left_best_hits, final_junctions); |
936 |
update_insertions_and_deletions(left_best_hits, final_insertions, final_deletions); |
937 |
|
938 |
print_sam_for_hits(rt, |
939 |
left_best_hits, |
940 |
grade, |
941 |
FRAG_LEFT, |
942 |
left_reads, |
943 |
bam_writer); |
944 |
} |
945 |
} |
946 |
else |
947 |
{ |
948 |
HitsForRead left_best_hits; |
949 |
HitsForRead right_best_hits; |
950 |
left_best_hits.insert_id = curr_left_obs_order; |
951 |
right_best_hits.insert_id = curr_right_obs_order; |
952 |
|
953 |
InsertAlignmentGrade grade; |
954 |
insert_best_alignments(curr_left_hit_group, |
955 |
curr_right_hit_group, |
956 |
grade, |
957 |
left_best_hits, |
958 |
right_best_hits); |
959 |
|
960 |
if (left_best_hits.hits.size() <= max_multihits && right_best_hits.hits.size() <= max_multihits) |
961 |
{ |
962 |
update_junctions(left_best_hits, final_junctions); |
963 |
update_junctions(right_best_hits, final_junctions); |
964 |
update_insertions_and_deletions(left_best_hits, final_insertions, final_deletions); |
965 |
update_insertions_and_deletions(right_best_hits, final_insertions, final_deletions); |
966 |
|
967 |
print_sam_for_hits(rt, |
968 |
left_best_hits, |
969 |
right_best_hits, |
970 |
grade, |
971 |
left_reads, |
972 |
right_reads, |
973 |
bam_writer); |
974 |
} |
975 |
} |
976 |
|
977 |
left_hs.next_read_hits(curr_left_hit_group); |
978 |
curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id); |
979 |
|
980 |
right_hs.next_read_hits(curr_right_hit_group); |
981 |
curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id); |
982 |
} |
983 |
|
984 |
} //while we still have unreported hits.. |
985 |
|
986 |
fprintf (stderr, "done.\n"); |
987 |
|
988 |
small_overhangs = 0; |
989 |
for (JunctionSet::iterator i = final_junctions.begin(); i != final_junctions.end();) |
990 |
{ |
991 |
if (i->second.supporting_hits == 0 || i->second.left_extent < 8 || i->second.right_extent < 8) |
992 |
{ |
993 |
final_junctions.erase(i++); |
994 |
} |
995 |
else |
996 |
{ |
997 |
++i; |
998 |
} |
999 |
} |
1000 |
|
1001 |
// if (small_overhangs > 0) |
1002 |
// fprintf(stderr, "Warning: %lu small overhang junctions!\n", small_overhangs); |
1003 |
|
1004 |
fprintf (stderr, "Printing junction BED track..."); |
1005 |
print_junctions(junctions_out, final_junctions, rt); |
1006 |
fprintf (stderr, "done\n"); |
1007 |
|
1008 |
fprintf (stderr, "Printing insertions..."); |
1009 |
print_insertions(insertions_out, final_insertions,rt); |
1010 |
fclose(insertions_out); |
1011 |
fprintf (stderr, "done\n"); |
1012 |
|
1013 |
fprintf (stderr, "Printing deletions..."); |
1014 |
print_deletions(deletions_out, final_deletions, rt); |
1015 |
fclose(deletions_out); |
1016 |
fprintf (stderr, "done\n"); |
1017 |
|
1018 |
fprintf(stderr, "Found %lu junctions from happy spliced reads\n", (long unsigned int)final_junctions.size()); |
1019 |
} |
1020 |
|
1021 |
void print_usage() |
1022 |
{ |
1023 |
fprintf(stderr, "Usage: tophat_reports <junctions.bed> <insertions.vcf> <deletions.vcf> <accepted_hits.sam> <left_map1,...,left_mapN> <left_reads.fq> [right_map1,...,right_mapN] [right_reads.fq]\n"); |
1024 |
|
1025 |
// fprintf(stderr, "Usage: tophat_reports <coverage.wig> <junctions.bed> <accepted_hits.sam> <map1.bwtout> [splice_map1.sbwtout]\n"); |
1026 |
} |
1027 |
|
1028 |
int main(int argc, char** argv) |
1029 |
{ |
1030 |
fprintf(stderr, "tophat_reports v%s (%s)\n", PACKAGE_VERSION, SVN_REVISION); |
1031 |
fprintf(stderr, "---------------------------------------\n"); |
1032 |
|
1033 |
reads_format = FASTQ; |
1034 |
|
1035 |
int parse_ret = parse_options(argc, argv, print_usage); |
1036 |
if (parse_ret) |
1037 |
return parse_ret; |
1038 |
|
1039 |
if(optind >= argc) |
1040 |
{ |
1041 |
print_usage(); |
1042 |
return 1; |
1043 |
} |
1044 |
|
1045 |
string junctions_file_name = argv[optind++]; |
1046 |
|
1047 |
if(optind >= argc) |
1048 |
{ |
1049 |
print_usage(); |
1050 |
return 1; |
1051 |
} |
1052 |
|
1053 |
string insertions_file_name = argv[optind++]; |
1054 |
if(optind >= argc) |
1055 |
{ |
1056 |
print_usage(); |
1057 |
return 1; |
1058 |
} |
1059 |
|
1060 |
string deletions_file_name = argv[optind++]; |
1061 |
if(optind >= argc) |
1062 |
{ |
1063 |
print_usage(); |
1064 |
return 1; |
1065 |
} |
1066 |
|
1067 |
string accepted_hits_file_name = argv[optind++]; |
1068 |
|
1069 |
if(optind >= argc) |
1070 |
{ |
1071 |
print_usage(); |
1072 |
return 1; |
1073 |
} |
1074 |
|
1075 |
string left_map_filename = argv[optind++]; |
1076 |
|
1077 |
if(optind >= argc) |
1078 |
{ |
1079 |
print_usage(); |
1080 |
return 1; |
1081 |
} |
1082 |
//FZPipe left_map_file; |
1083 |
//string unbamcmd=getBam2SamCmd(left_map_filename); |
1084 |
//left_map_file.openRead(left_map_filename, unbamcmd); |
1085 |
string left_reads_filename = argv[optind++]; |
1086 |
string unzcmd=getUnpackCmd(left_reads_filename, false); |
1087 |
|
1088 |
string right_map_filename; |
1089 |
string right_reads_filename; |
1090 |
FZPipe right_reads_file; |
1091 |
|
1092 |
if (optind < argc) |
1093 |
{ |
1094 |
right_map_filename = argv[optind++]; |
1095 |
if(optind >= argc) { |
1096 |
print_usage(); |
1097 |
return 1; |
1098 |
} |
1099 |
//right_map_file.openRead(right_map_filename, unbamcmd); |
1100 |
//if (optind<argc) { |
1101 |
right_reads_filename=argv[optind++]; |
1102 |
right_reads_file.openRead(right_reads_filename,unzcmd); |
1103 |
// } |
1104 |
} |
1105 |
FILE* junctions_file = fopen(junctions_file_name.c_str(), "w"); |
1106 |
if (junctions_file == NULL) |
1107 |
{ |
1108 |
fprintf(stderr, "Error: cannot open BED file %s for writing\n", |
1109 |
junctions_file_name.c_str()); |
1110 |
exit(1); |
1111 |
} |
1112 |
|
1113 |
FILE* insertions_file = fopen(insertions_file_name.c_str(), "w"); |
1114 |
if (insertions_file == NULL) |
1115 |
{ |
1116 |
fprintf(stderr, "Error: cannot open VCF file %s for writing\n", |
1117 |
insertions_file_name.c_str()); |
1118 |
exit(1); |
1119 |
} |
1120 |
|
1121 |
|
1122 |
FILE* deletions_file = fopen(deletions_file_name.c_str(), "w"); |
1123 |
if (deletions_file == NULL) |
1124 |
{ |
1125 |
fprintf(stderr, "Error: cannot open VCF file %s for writing\n", |
1126 |
deletions_file_name.c_str()); |
1127 |
exit(1); |
1128 |
} |
1129 |
bool uncompressed_bam=(accepted_hits_file_name=="-"); |
1130 |
GBamWriter bam_writer(accepted_hits_file_name.c_str(), sam_header.c_str(), uncompressed_bam); |
1131 |
FZPipe left_reads_file(left_reads_filename, unzcmd); |
1132 |
if (left_reads_file.file==NULL) |
1133 |
{ |
1134 |
fprintf(stderr, "Error: cannot open reads file %s for reading\n", |
1135 |
left_reads_filename.c_str()); |
1136 |
exit(1); |
1137 |
} |
1138 |
|
1139 |
driver(bam_writer, left_map_filename, |
1140 |
left_reads_file.file, |
1141 |
right_map_filename, |
1142 |
right_reads_file.file, |
1143 |
junctions_file, |
1144 |
insertions_file, |
1145 |
deletions_file); |
1146 |
return 0; |
1147 |
} |