26 |
|
#include <seqan/sequence.h> |
27 |
|
#include <seqan/find.h> |
28 |
|
#include <seqan/file.h> |
29 |
+ |
#include <seqan/modifier.h> |
30 |
|
#include <getopt.h> |
31 |
|
|
32 |
|
#include "common.h" |
35 |
|
#include "junctions.h" |
36 |
|
#include "insertions.h" |
37 |
|
#include "deletions.h" |
38 |
+ |
#include "fusions.h" |
39 |
+ |
|
40 |
|
|
41 |
|
using namespace std; |
42 |
|
using namespace seqan; |
101 |
|
|
102 |
|
template<typename TStr> |
103 |
|
void print_splice(const Junction& junction, |
104 |
< |
int read_len, |
105 |
< |
const string& tag, |
106 |
< |
TStr& ref_str, |
107 |
< |
const string& ref_name, |
108 |
< |
ostream& splice_db) |
104 |
> |
int read_len, |
105 |
> |
const string& tag, |
106 |
> |
TStr& ref_str, |
107 |
> |
const string& ref_name, |
108 |
> |
ostream& splice_db) |
109 |
|
{ |
110 |
|
// daehwan - this is tentative, let's think about this more :) |
111 |
|
// int half_splice_len = read_len - min_anchor_len; |
112 |
|
int half_splice_len = read_len; |
110 |
– |
|
113 |
|
size_t left_start, right_start; |
114 |
|
size_t left_end, right_end; |
115 |
|
|
136 |
|
} |
137 |
|
} |
138 |
|
|
139 |
+ |
template<typename TStr> |
140 |
+ |
void print_fusion(const Fusion& fusion, |
141 |
+ |
int read_len, |
142 |
+ |
TStr& left_ref_str, |
143 |
+ |
TStr& right_ref_str, |
144 |
+ |
const char* left_ref_name, |
145 |
+ |
const char* right_ref_name, |
146 |
+ |
ostream& fusion_db) |
147 |
+ |
{ |
148 |
+ |
int half_splice_len = read_len - min_anchor_len; |
149 |
+ |
|
150 |
+ |
size_t left_start, right_start; |
151 |
+ |
size_t left_end, right_end; |
152 |
+ |
|
153 |
+ |
if (fusion.left >= 0 && fusion.left < length(left_ref_str) && |
154 |
+ |
fusion.right >= 0 && fusion.right < length(right_ref_str)) |
155 |
+ |
{ |
156 |
+ |
if (fusion.dir == FUSION_FF || fusion.dir == FUSION_FR) |
157 |
+ |
{ |
158 |
+ |
left_start = fusion.left + 1 >= half_splice_len ? fusion.left - half_splice_len + 1 : 0; |
159 |
+ |
left_end = left_start + half_splice_len; |
160 |
+ |
} |
161 |
+ |
else |
162 |
+ |
{ |
163 |
+ |
left_start = fusion.left; |
164 |
+ |
left_end = left_start + half_splice_len < length(left_ref_str) ? left_start + half_splice_len : length(left_ref_str); |
165 |
+ |
} |
166 |
+ |
|
167 |
+ |
if (fusion.dir == FUSION_FF || fusion.dir == FUSION_RF) |
168 |
+ |
{ |
169 |
+ |
right_start = fusion.right; |
170 |
+ |
right_end = right_start + half_splice_len < length(right_ref_str) ? right_start + half_splice_len : length(right_ref_str); |
171 |
+ |
} |
172 |
+ |
else |
173 |
+ |
{ |
174 |
+ |
right_end = fusion.right + 1; |
175 |
+ |
right_start = right_end >= half_splice_len ? right_end - half_splice_len : 0; |
176 |
+ |
} |
177 |
+ |
|
178 |
+ |
seqan::Dna5String left_splice = infix(left_ref_str, left_start, left_end); |
179 |
+ |
seqan::Dna5String right_splice = infix(right_ref_str, right_start, right_end); |
180 |
+ |
|
181 |
+ |
if (fusion.dir == FUSION_RF || fusion.dir == FUSION_RR) |
182 |
+ |
{ |
183 |
+ |
seqan::convertInPlace(left_splice, seqan::FunctorComplement<Dna>()); |
184 |
+ |
seqan::reverseInPlace(left_splice); |
185 |
+ |
|
186 |
+ |
left_start = left_end - 1; |
187 |
+ |
} |
188 |
+ |
|
189 |
+ |
if (fusion.dir == FUSION_FR || fusion.dir == FUSION_RR) |
190 |
+ |
{ |
191 |
+ |
seqan::convertInPlace(right_splice, seqan::FunctorComplement<Dna>()); |
192 |
+ |
seqan::reverseInPlace(right_splice); |
193 |
+ |
|
194 |
+ |
right_end = right_start - 1; |
195 |
+ |
} |
196 |
+ |
|
197 |
+ |
const char* dir = "ff"; |
198 |
+ |
if (fusion.dir == FUSION_FR) |
199 |
+ |
dir = "fr"; |
200 |
+ |
else if (fusion.dir == FUSION_RF) |
201 |
+ |
dir = "rf"; |
202 |
+ |
else if (fusion.dir == FUSION_RR) |
203 |
+ |
dir = "rr"; |
204 |
+ |
|
205 |
+ |
fusion_db << ">" << left_ref_name << "-" << right_ref_name << "|" |
206 |
+ |
<< left_start << "|" |
207 |
+ |
<< fusion.left << "-" << fusion.right << "|" |
208 |
+ |
<< right_end << "|fus|" << dir << endl; |
209 |
+ |
|
210 |
+ |
fusion_db << left_splice << right_splice << endl; |
211 |
+ |
} |
212 |
+ |
} |
213 |
+ |
|
214 |
|
|
215 |
|
/** |
216 |
|
* Parse an int out of optarg and enforce that it be at least 'lower'; |
256 |
|
// return 0; |
257 |
|
//} |
258 |
|
|
259 |
+ |
void get_seqs(istream& ref_stream, |
260 |
+ |
RefSequenceTable& rt, |
261 |
+ |
bool keep_seqs = true) |
262 |
+ |
{ |
263 |
+ |
while(ref_stream.good() && |
264 |
+ |
!ref_stream.eof()) |
265 |
+ |
{ |
266 |
+ |
RefSequenceTable::Sequence* ref_str = new RefSequenceTable::Sequence(); |
267 |
+ |
string name; |
268 |
+ |
readMeta(ref_stream, name, Fasta()); |
269 |
+ |
string::size_type space_pos = name.find_first_of(" \t\r"); |
270 |
+ |
if (space_pos != string::npos) |
271 |
+ |
{ |
272 |
+ |
name.resize(space_pos); |
273 |
+ |
} |
274 |
+ |
fprintf(stderr, "\tLoading %s...", name.c_str()); |
275 |
+ |
seqan::read(ref_stream, *ref_str, Fasta()); |
276 |
+ |
fprintf(stderr, "done\n"); |
277 |
+ |
rt.get_id(name, keep_seqs ? ref_str : NULL, 0); |
278 |
+ |
if (!keep_seqs) |
279 |
+ |
delete ref_str; |
280 |
+ |
} |
281 |
+ |
} |
282 |
+ |
|
283 |
|
void driver(const vector<FILE*>& splice_coords_files, |
284 |
< |
const vector<FILE*>& insertion_coords_files, |
285 |
< |
const vector<FILE*>& deletion_coords_files, |
286 |
< |
ifstream& ref_stream) |
284 |
> |
const vector<FILE*>& insertion_coords_files, |
285 |
> |
const vector<FILE*>& deletion_coords_files, |
286 |
> |
const vector<FILE*>& fusion_coords_files, |
287 |
> |
ifstream& ref_stream) |
288 |
|
{ |
289 |
|
char splice_buf[2048]; |
290 |
< |
RefSequenceTable rt(true); |
290 |
> |
RefSequenceTable rt(sam_header, true); |
291 |
> |
get_seqs(ref_stream, rt, true); |
292 |
> |
|
293 |
|
JunctionSet junctions; |
294 |
|
for (size_t i = 0; i < splice_coords_files.size(); ++i) |
295 |
|
{ |
420 |
|
} |
421 |
|
} |
422 |
|
|
423 |
+ |
std::set<Fusion> fusions; |
424 |
+ |
for(size_t i=0; i < fusion_coords_files.size(); ++i){ |
425 |
+ |
FILE* fusion_coords = fusion_coords_files[i]; |
426 |
+ |
if(!fusion_coords){ |
427 |
+ |
continue; |
428 |
+ |
} |
429 |
+ |
while(fgets(splice_buf, 2048, fusion_coords)){ |
430 |
+ |
char* nl = strrchr(splice_buf, '\n'); |
431 |
+ |
char* buf = splice_buf; |
432 |
+ |
if (nl) *nl = 0; |
433 |
+ |
|
434 |
+ |
char* ref_name1 = strsep((char**)&buf, "\t"); |
435 |
+ |
char* scan_left_coord = strsep((char**)&buf, "\t"); |
436 |
+ |
char* ref_name2 = strsep((char**)&buf, "\t"); |
437 |
+ |
char* scan_right_coord = strsep((char**)&buf, "\t"); |
438 |
+ |
char* scan_dir = strsep((char**)&buf, "\t"); |
439 |
|
|
440 |
< |
typedef RefSequenceTable::Sequence Reference; |
441 |
< |
|
442 |
< |
while(ref_stream.good() && |
443 |
< |
!ref_stream.eof()) |
444 |
< |
{ |
445 |
< |
Reference ref_str; |
446 |
< |
string name; |
447 |
< |
|
448 |
< |
readMeta(ref_stream, name, Fasta()); |
449 |
< |
string::size_type space_pos = name.find_first_of(" \t\r"); |
450 |
< |
if (space_pos != string::npos) |
451 |
< |
{ |
452 |
< |
name.resize(space_pos); |
453 |
< |
} |
454 |
< |
|
455 |
< |
read(ref_stream, ref_str, Fasta()); |
456 |
< |
|
457 |
< |
uint32_t refid = rt.get_id(name, NULL, 0); |
458 |
< |
Junction dummy_left(refid, 0, 0, true); |
339 |
< |
Junction dummy_right(refid, VMAXINT32, VMAXINT32, true); |
340 |
< |
|
341 |
< |
pair<JunctionSet::iterator, JunctionSet::iterator> r; |
342 |
< |
r.first = junctions.lower_bound(dummy_left); |
343 |
< |
r.second = junctions.upper_bound(dummy_right); |
344 |
< |
|
345 |
< |
JunctionSet::iterator itr = r.first; |
346 |
< |
|
347 |
< |
while(itr != r.second && itr != junctions.end()) |
348 |
< |
{ |
349 |
< |
print_splice(itr->first, read_length, itr->first.antisense ? "GTAG|rev" : "GTAG|fwd", ref_str, name, cout); |
350 |
< |
++itr; |
440 |
> |
if (!ref_name1 || !scan_left_coord || !ref_name2 || !scan_right_coord || !scan_dir) |
441 |
> |
{ |
442 |
> |
fprintf(stderr,"Error: malformed insertion coordinate record\n"); |
443 |
> |
exit(1); |
444 |
> |
} |
445 |
> |
|
446 |
> |
uint32_t ref_id1 = rt.get_id(ref_name1, NULL, 0); |
447 |
> |
uint32_t ref_id2 = rt.get_id(ref_name2, NULL, 0); |
448 |
> |
uint32_t left_coord = atoi(scan_left_coord); |
449 |
> |
uint32_t right_coord = atoi(scan_right_coord); |
450 |
> |
uint32_t dir = FUSION_FF; |
451 |
> |
if (strcmp(scan_dir, "fr") == 0) |
452 |
> |
dir = FUSION_FR; |
453 |
> |
else if(strcmp(scan_dir, "rf") == 0) |
454 |
> |
dir = FUSION_RF; |
455 |
> |
else if(strcmp(scan_dir, "rr") == 0) |
456 |
> |
dir = FUSION_RR; |
457 |
> |
|
458 |
> |
fusions.insert(Fusion(ref_id1, ref_id2, left_coord, right_coord, dir)); |
459 |
|
} |
460 |
|
} |
461 |
|
|
354 |
– |
|
355 |
– |
ref_stream.clear(); |
356 |
– |
ref_stream.seekg(0, ios::beg); |
357 |
– |
|
358 |
– |
|
359 |
– |
while(ref_stream.good() && |
360 |
– |
!ref_stream.eof()) |
462 |
|
{ |
463 |
< |
Reference ref_str; |
464 |
< |
string name; |
465 |
< |
|
466 |
< |
readMeta(ref_stream, name, Fasta()); |
467 |
< |
string::size_type space_pos = name.find_first_of(" \t\r"); |
468 |
< |
if (space_pos != string::npos) |
469 |
< |
{ |
470 |
< |
name.resize(space_pos); |
370 |
< |
} |
371 |
< |
|
372 |
< |
read(ref_stream, ref_str, Fasta()); |
373 |
< |
|
374 |
< |
uint32_t refid = rt.get_id(name, NULL,0); |
375 |
< |
Deletion dummy_left(refid, 0, 0, true); |
376 |
< |
Deletion dummy_right(refid, VMAXINT32, VMAXINT32, true); |
377 |
< |
|
378 |
< |
pair<std::set<Deletion>::iterator, std::set<Deletion>::iterator> r; |
379 |
< |
r.first = deletions.lower_bound(dummy_left); |
380 |
< |
r.second = deletions.upper_bound(dummy_right); |
381 |
< |
|
382 |
< |
std::set<Deletion>::iterator itr = r.first; |
383 |
< |
|
384 |
< |
while(itr != r.second && itr != deletions.end()) |
385 |
< |
{ |
386 |
< |
print_splice((Junction)*itr, read_length, itr->antisense ? "del|rev" : "del|fwd", ref_str, name, cout); |
387 |
< |
++itr; |
388 |
< |
} |
463 |
> |
JunctionSet::iterator itr = junctions.begin(); |
464 |
> |
while(itr != junctions.end()) |
465 |
> |
{ |
466 |
> |
RefSequenceTable::Sequence* ref_str = rt.get_seq(itr->first.refid); |
467 |
> |
const char* name = rt.get_name(itr->first.refid); |
468 |
> |
print_splice(itr->first, read_length, itr->first.antisense ? "GTAG|rev" : "GTAG|fwd", *ref_str, name, cout); |
469 |
> |
++itr; |
470 |
> |
} |
471 |
|
} |
472 |
|
|
391 |
– |
ref_stream.clear(); |
392 |
– |
ref_stream.seekg(0, ios::beg); |
393 |
– |
|
394 |
– |
|
395 |
– |
|
396 |
– |
while(ref_stream.good() && |
397 |
– |
!ref_stream.eof()) |
473 |
|
{ |
474 |
< |
Reference ref_str; |
475 |
< |
string name; |
476 |
< |
|
477 |
< |
readMeta(ref_stream, name, Fasta()); |
478 |
< |
string::size_type space_pos = name.find_first_of(" \t\r"); |
479 |
< |
if (space_pos != string::npos) |
480 |
< |
{ |
481 |
< |
name.resize(space_pos); |
482 |
< |
} |
408 |
< |
|
409 |
< |
read(ref_stream, ref_str, Fasta()); |
410 |
< |
|
411 |
< |
uint32_t refid = rt.get_id(name, NULL,0); |
412 |
< |
Insertion dummy_left(refid, 0, ""); |
413 |
< |
Insertion dummy_right(refid, VMAXINT32, ""); |
414 |
< |
|
415 |
< |
std::set<Insertion>::iterator itr = insertions.lower_bound(dummy_left); |
416 |
< |
std::set<Insertion>::iterator upper = insertions.upper_bound(dummy_right); |
474 |
> |
std::set<Deletion>::iterator itr = deletions.begin(); |
475 |
> |
while(itr != deletions.end()) |
476 |
> |
{ |
477 |
> |
RefSequenceTable::Sequence* ref_str = rt.get_seq(itr->refid); |
478 |
> |
const char* name = rt.get_name(itr->refid); |
479 |
> |
print_splice((Junction)*itr, read_length, itr->antisense ? "del|rev" : "del|fwd", *ref_str, name, cout); |
480 |
> |
++itr; |
481 |
> |
} |
482 |
> |
} |
483 |
|
|
484 |
< |
while(itr != upper && itr != insertions.end()){ |
485 |
< |
print_insertion(*itr, read_length, ref_str, name, cout); |
486 |
< |
++itr; |
487 |
< |
} |
484 |
> |
{ |
485 |
> |
std::set<Insertion>::iterator itr = insertions.begin(); |
486 |
> |
while(itr != insertions.end()){ |
487 |
> |
RefSequenceTable::Sequence* ref_str = rt.get_seq(itr->refid); |
488 |
> |
const char* name = rt.get_name(itr->refid); |
489 |
> |
print_insertion(*itr, read_length, *ref_str, name, cout); |
490 |
> |
++itr; |
491 |
> |
} |
492 |
|
} |
493 |
|
|
494 |
+ |
{ |
495 |
+ |
std::set<Fusion>::iterator itr = fusions.begin(); |
496 |
+ |
while(itr != fusions.end()){ |
497 |
+ |
RefSequenceTable::Sequence* left_ref_str = rt.get_seq(itr->refid1); |
498 |
+ |
RefSequenceTable::Sequence* right_ref_str = rt.get_seq(itr->refid2); |
499 |
+ |
const char* left_ref_name = rt.get_name(itr->refid1); |
500 |
+ |
const char* right_ref_name = rt.get_name(itr->refid2); |
501 |
+ |
print_fusion(*itr, read_length, *left_ref_str, *right_ref_str, left_ref_name, right_ref_name, cout); |
502 |
+ |
++itr; |
503 |
+ |
} |
504 |
+ |
} |
505 |
|
} |
506 |
|
|
507 |
|
int main(int argc, char** argv) |
608 |
|
} |
609 |
|
|
610 |
|
|
611 |
+ |
/* |
612 |
+ |
*/ |
613 |
+ |
string fusion_coords_file_list = argv[optind++]; |
614 |
+ |
vector<string> fusion_coords_file_names; |
615 |
+ |
vector<FILE*> fusion_coords_files; |
616 |
+ |
tokenize(fusion_coords_file_list, ",", fusion_coords_file_names); |
617 |
+ |
for(size_t s = 0; s < fusion_coords_file_names.size(); ++s) |
618 |
+ |
{ |
619 |
+ |
FILE* fusion_coords_file = fopen(fusion_coords_file_names[s].c_str(),"r"); |
620 |
+ |
if(!fusion_coords_file) |
621 |
+ |
{ |
622 |
+ |
fprintf(stderr, "Warning: cannot open %s for reading\n", |
623 |
+ |
fusion_coords_file_names[s].c_str()); |
624 |
+ |
continue; |
625 |
+ |
} |
626 |
+ |
fusion_coords_files.push_back(fusion_coords_file); |
627 |
+ |
} |
628 |
+ |
if(optind >= argc) |
629 |
+ |
{ |
630 |
+ |
print_usage(); |
631 |
+ |
return 1; |
632 |
+ |
} |
633 |
+ |
|
634 |
|
|
635 |
|
string ref_file_name = argv[optind++]; |
636 |
|
ifstream ref_stream(ref_file_name.c_str()); |
642 |
|
exit(1); |
643 |
|
} |
644 |
|
|
645 |
< |
driver(coords_files, insertion_coords_files, deletion_coords_files, ref_stream); |
645 |
> |
driver(coords_files, insertion_coords_files, deletion_coords_files, fusion_coords_files, ref_stream); |
646 |
|
return 0; |
647 |
|
} |