14 |
|
#include <iostream> |
15 |
|
#include <fstream> |
16 |
|
#include <getopt.h> |
17 |
+ |
#include <seqan/sequence.h> |
18 |
+ |
#include <seqan/file.h> |
19 |
+ |
|
20 |
|
#include "bwt_map.h" |
21 |
|
#include "inserts.h" |
22 |
|
#include "closures.h" |
23 |
|
#include "reads.h" |
24 |
|
#include "tokenize.h" |
25 |
< |
#include <seqan/sequence.h> |
23 |
< |
#include <seqan/file.h> |
25 |
> |
#include "fusions.h" |
26 |
|
|
27 |
|
using namespace seqan; |
28 |
|
using namespace std; |
86 |
|
} |
87 |
|
} |
88 |
|
|
89 |
+ |
void find_fusion_closure(HitsForRead& left_hits, |
90 |
+ |
HitsForRead& right_hits, |
91 |
+ |
std::set<Fusion>& fusions) |
92 |
+ |
{ |
93 |
+ |
if (left_hits.hits.size() > 3 || right_hits.hits.size() > 3) |
94 |
+ |
return; |
95 |
+ |
|
96 |
+ |
for (size_t left_index = 0; left_index < left_hits.hits.size(); ++left_index) |
97 |
+ |
{ |
98 |
+ |
for (size_t right_index = 0; right_index < right_hits.hits.size(); ++right_index) |
99 |
+ |
{ |
100 |
+ |
BowtieHit* leftHit = &left_hits.hits[left_index]; |
101 |
+ |
BowtieHit* rightHit = &right_hits.hits[right_index]; |
102 |
+ |
|
103 |
+ |
uint32_t dir = FUSION_FF; |
104 |
+ |
if (leftHit->antisense_align() != rightHit->antisense_align()) |
105 |
+ |
dir = FUSION_FF; |
106 |
+ |
else if (!leftHit->antisense_align()) |
107 |
+ |
dir = FUSION_FR; |
108 |
+ |
else |
109 |
+ |
dir = FUSION_RF; |
110 |
+ |
|
111 |
+ |
if (dir == FUSION_FF && leftHit->antisense_align()) |
112 |
+ |
{ |
113 |
+ |
BowtieHit * tmp = leftHit; |
114 |
+ |
leftHit = rightHit; |
115 |
+ |
rightHit = tmp; |
116 |
+ |
} |
117 |
+ |
|
118 |
+ |
uint32_t left, right; |
119 |
+ |
if (dir == FUSION_FF || dir == FUSION_FR) |
120 |
+ |
left = leftHit->right() - 4; |
121 |
+ |
else |
122 |
+ |
left = leftHit->left() + 4; |
123 |
+ |
|
124 |
+ |
if (dir == FUSION_FF || dir == FUSION_RF) |
125 |
+ |
right = rightHit->left() + 4; |
126 |
+ |
else |
127 |
+ |
right = rightHit->right() - 4; |
128 |
+ |
|
129 |
+ |
if (leftHit->ref_id() == rightHit->ref_id() && dir == FUSION_FF) |
130 |
+ |
{ |
131 |
+ |
int dist = 0; |
132 |
+ |
dist = rightHit->left() - leftHit->right(); |
133 |
+ |
|
134 |
+ |
if (dist >= 0 && dist <= (int)fusion_min_dist) |
135 |
+ |
continue; |
136 |
+ |
} |
137 |
+ |
|
138 |
+ |
uint32_t ref_id1 = leftHit->ref_id(); |
139 |
+ |
uint32_t ref_id2 = rightHit->ref_id(); |
140 |
+ |
|
141 |
+ |
if (dir == FUSION_FR || dir == FUSION_RF) |
142 |
+ |
{ |
143 |
+ |
if ((ref_id2 < ref_id1) || (ref_id1 == ref_id2 && left > right)) |
144 |
+ |
{ |
145 |
+ |
uint32_t temp = ref_id1; |
146 |
+ |
ref_id1 = ref_id2; |
147 |
+ |
ref_id2 = temp; |
148 |
+ |
|
149 |
+ |
temp = left; |
150 |
+ |
left = right; |
151 |
+ |
right = temp; |
152 |
+ |
} |
153 |
+ |
} |
154 |
+ |
|
155 |
+ |
fusions.insert(Fusion(ref_id1, ref_id2, left, right, dir)); |
156 |
+ |
} |
157 |
+ |
} |
158 |
+ |
} |
159 |
+ |
|
160 |
|
bool prefer_shorter_pairs = true; |
161 |
|
|
162 |
|
typedef pair<InsertAlignmentGrade, vector<InsertAlignment> > BestPairingHits; |
513 |
|
void closure_driver(vector<FZPipe>& map1, |
514 |
|
vector<FZPipe>& map2, |
515 |
|
ifstream& ref_stream, |
516 |
< |
FILE* juncs_file) |
516 |
> |
FILE* juncs_file, |
517 |
> |
FILE* fusions_out) |
518 |
|
{ |
519 |
|
typedef RefSequenceTable::Sequence Reference; |
520 |
|
|
521 |
|
ReadTable it; |
522 |
|
RefSequenceTable rt(true); |
523 |
< |
|
524 |
< |
BowtieHitFactory hit_factory(it,rt); |
523 |
> |
|
524 |
> |
BowtieHitFactory hit_factory(it, rt); |
525 |
> |
|
526 |
> |
std::set<Fusion> fusions; |
527 |
|
|
528 |
|
fprintf (stderr, "Finding near-covered motifs..."); |
529 |
|
CoverageMapVisitor cov_map_visitor(ref_stream, rt); |
567 |
|
|
568 |
|
while (curr_left_obs_order == curr_right_obs_order && |
569 |
|
curr_left_obs_order != VMAXINT32 && curr_right_obs_order != VMAXINT32) |
570 |
< |
{ |
571 |
< |
if (coverage_attempts++ % 1000 == 0) |
572 |
< |
fprintf (stderr, "Adding covered motifs from pair %d\n", coverage_attempts); |
570 |
> |
{ |
571 |
> |
if (num == 0) |
572 |
> |
find_fusion_closure(curr_left_hit_group, curr_right_hit_group, fusions); |
573 |
> |
|
574 |
> |
if (coverage_attempts++ % 10000 == 0) |
575 |
> |
fprintf (stderr, "Adding covered motifs from pair %d\n", coverage_attempts); |
576 |
> |
|
577 |
|
visit_best_pairing(curr_left_hit_group, curr_right_hit_group, cov_map_visitor); |
578 |
|
|
579 |
|
left_hs.next_read_hits(curr_left_hit_group); |
636 |
|
while (curr_left_obs_order == curr_right_obs_order && |
637 |
|
curr_left_obs_order != VMAXINT32 && curr_right_obs_order != VMAXINT32) |
638 |
|
{ |
639 |
< |
if (closure_attempts++ % 1000 == 0) |
640 |
< |
fprintf (stderr, "Trying to close pair %d\n", closure_attempts); |
639 |
> |
if (closure_attempts++ % 10000 == 0) |
640 |
> |
fprintf (stderr, "Trying to close pair %d\n", closure_attempts); |
641 |
> |
|
642 |
|
visit_best_pairing(curr_left_hit_group, curr_right_hit_group, junc_map_visitor); |
643 |
|
left_hs.next_read_hits(curr_left_hit_group); |
644 |
|
curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id); |
691 |
|
fprintf(stderr, "Successfully closed %d pairs\n", closed); |
692 |
|
|
693 |
|
fprintf(stderr, "Found %d total possible splices\n", num_potential_splices); |
694 |
+ |
|
695 |
+ |
// daehwan |
696 |
+ |
#if 0 |
697 |
+ |
fprintf (stderr, "Reporting potential fusions...\n"); |
698 |
+ |
if(fusions_out){ |
699 |
+ |
for(std::set<Fusion>::iterator itr = fusions.begin(); itr != fusions.end(); ++itr){ |
700 |
+ |
const char* ref_name1 = rt.get_name(itr->refid1); |
701 |
+ |
const char* ref_name2 = rt.get_name(itr->refid2); |
702 |
+ |
|
703 |
+ |
const char* dir = ""; |
704 |
+ |
if (itr->dir == FUSION_FR) |
705 |
+ |
dir = "fr"; |
706 |
+ |
else if(itr->dir == FUSION_RF) |
707 |
+ |
dir = "rf"; |
708 |
+ |
else |
709 |
+ |
dir = "ff"; |
710 |
+ |
|
711 |
+ |
fprintf(fusions_out, |
712 |
+ |
"%s\t%d\t%s\t%d\t%s\n", |
713 |
+ |
ref_name1, |
714 |
+ |
itr->left, |
715 |
+ |
ref_name2, |
716 |
+ |
itr->right, |
717 |
+ |
dir); |
718 |
+ |
} |
719 |
+ |
fclose(fusions_out); |
720 |
+ |
}else{ |
721 |
+ |
fprintf(stderr, "Failed to open fusions file for writing\n"); |
722 |
+ |
} |
723 |
+ |
#endif |
724 |
|
} |
725 |
|
|
726 |
|
void print_usage() |
750 |
|
print_usage(); |
751 |
|
return 1; |
752 |
|
} |
753 |
+ |
|
754 |
+ |
string fusions_file_name = argv[optind++]; |
755 |
|
|
756 |
+ |
if(optind >= argc) |
757 |
+ |
{ |
758 |
+ |
print_usage(); |
759 |
+ |
return 1; |
760 |
+ |
} |
761 |
+ |
|
762 |
|
string ref_fasta = argv[optind++]; |
763 |
< |
|
763 |
> |
|
764 |
|
if(optind >= argc) |
765 |
|
{ |
766 |
|
print_usage(); |
816 |
|
junctions_file_name.c_str()); |
817 |
|
exit(1); |
818 |
|
} |
819 |
< |
|
819 |
> |
|
820 |
> |
FILE* fusion_db = fopen(fusions_file_name.c_str(), "w"); |
821 |
> |
if (splice_db == NULL) |
822 |
> |
{ |
823 |
> |
fprintf(stderr, "Error: cannot open fusions file %s for writing\n", |
824 |
> |
fusions_file_name.c_str()); |
825 |
> |
exit(1); |
826 |
> |
} |
827 |
> |
|
828 |
|
closure_driver(left_files, |
829 |
|
right_files, |
830 |
|
ref_stream, |
831 |
< |
splice_db); |
831 |
> |
splice_db, |
832 |
> |
fusion_db); |
833 |
|
|
834 |
|
return 0; |
835 |
|
} |