73 |
|
} |
74 |
|
} |
75 |
|
|
76 |
+ |
struct cmp_read_alignment |
77 |
+ |
{ |
78 |
+ |
cmp_read_alignment(const JunctionSet& gtf_juncs) : |
79 |
+ |
_gtf_juncs(>f_juncs) {} |
80 |
+ |
|
81 |
+ |
bool operator() (const BowtieHit& l, const BowtieHit& r) const |
82 |
+ |
{ |
83 |
+ |
FragmentAlignmentGrade gl(l, *_gtf_juncs); |
84 |
+ |
FragmentAlignmentGrade gr(r, *_gtf_juncs); |
85 |
+ |
bool better = gr < gl; |
86 |
+ |
bool worse = gl < gr; |
87 |
+ |
|
88 |
+ |
if (better && !worse) |
89 |
+ |
return true; |
90 |
+ |
else |
91 |
+ |
return false; |
92 |
+ |
} |
93 |
+ |
|
94 |
+ |
const JunctionSet* _gtf_juncs; |
95 |
+ |
}; |
96 |
+ |
|
97 |
|
void read_best_alignments(const HitsForRead& hits_for_read, |
98 |
< |
FragmentAlignmentGrade& best_grade, |
99 |
< |
HitsForRead& best_hits, |
100 |
< |
const JunctionSet& gtf_junctions) |
98 |
> |
FragmentAlignmentGrade& best_grade, |
99 |
> |
HitsForRead& best_hits, |
100 |
> |
const JunctionSet& gtf_junctions) |
101 |
|
{ |
102 |
|
const vector<BowtieHit>& hits = hits_for_read.hits; |
103 |
|
for (size_t i = 0; i < hits.size(); ++i) |
104 |
|
{ |
105 |
< |
if (hits[i].edit_dist()>max_read_mismatches) continue; |
106 |
< |
|
107 |
< |
FragmentAlignmentGrade g(hits[i], gtf_junctions); |
108 |
< |
|
109 |
< |
// Is the new status better than the current best one? |
110 |
< |
if (best_grade < g) |
111 |
< |
{ |
112 |
< |
best_hits.hits.clear(); |
113 |
< |
best_grade = g; |
114 |
< |
best_hits.hits.push_back(hits[i]); |
115 |
< |
} |
116 |
< |
else if (!(g < best_grade)) // is it just as good? |
117 |
< |
{ |
118 |
< |
best_grade.num_alignments++; |
119 |
< |
best_hits.hits.push_back(hits[i]); |
120 |
< |
} |
105 |
> |
if (hits[i].edit_dist() > max_read_mismatches) |
106 |
> |
continue; |
107 |
> |
|
108 |
> |
if (report_secondary_alignments) |
109 |
> |
{ |
110 |
> |
best_hits.hits.push_back(hits[i]); |
111 |
> |
} |
112 |
> |
else |
113 |
> |
{ |
114 |
> |
FragmentAlignmentGrade g(hits[i], gtf_junctions); |
115 |
> |
|
116 |
> |
// Is the new status better than the current best one? |
117 |
> |
if (best_grade < g) |
118 |
> |
{ |
119 |
> |
best_hits.hits.clear(); |
120 |
> |
best_grade = g; |
121 |
> |
best_hits.hits.push_back(hits[i]); |
122 |
> |
} |
123 |
> |
else if (!(g < best_grade)) // is it just as good? |
124 |
> |
{ |
125 |
> |
best_grade.num_alignments++; |
126 |
> |
best_hits.hits.push_back(hits[i]); |
127 |
> |
} |
128 |
> |
} |
129 |
> |
} |
130 |
> |
|
131 |
> |
if (report_secondary_alignments && best_hits.hits.size() > 0) |
132 |
> |
{ |
133 |
> |
cmp_read_alignment cmp(gtf_junctions); |
134 |
> |
sort(best_hits.hits.begin(), best_hits.hits.end(), cmp); |
135 |
> |
best_grade = FragmentAlignmentGrade(best_hits.hits[0], gtf_junctions); |
136 |
> |
best_grade.num_alignments = best_hits.hits.size(); |
137 |
|
} |
138 |
|
} |
139 |
|
|
140 |
< |
void pair_best_alignments(const HitsForRead& left_hits, |
104 |
< |
const HitsForRead& right_hits, |
105 |
< |
InsertAlignmentGrade& best_grade, |
106 |
< |
HitsForRead& left_best_hits, |
107 |
< |
HitsForRead& right_best_hits) |
140 |
> |
void set_insert_alignment_grade(const BowtieHit& lh, const BowtieHit& rh, InsertAlignmentGrade& grade) |
141 |
|
{ |
142 |
|
// max mate inner distance (genomic) |
143 |
|
int min_mate_inner_dist = inner_dist_mean - inner_dist_std_dev; |
144 |
|
if (max_mate_inner_dist == -1) |
145 |
|
max_mate_inner_dist = inner_dist_mean + inner_dist_std_dev; |
146 |
|
|
147 |
+ |
bool fusion = false; |
148 |
+ |
if (fusion_search) |
149 |
+ |
{ |
150 |
+ |
bool left_fusion = lh.fusion_opcode() != FUSION_NOTHING; |
151 |
+ |
bool right_fusion = rh.fusion_opcode() != FUSION_NOTHING; |
152 |
+ |
if (left_fusion && right_fusion) |
153 |
+ |
return; |
154 |
+ |
|
155 |
+ |
fusion = left_fusion || right_fusion; |
156 |
+ |
if (!fusion && lh.ref_id() != rh.ref_id()) |
157 |
+ |
fusion = true; |
158 |
+ |
|
159 |
+ |
if (!fusion && lh.ref_id() == rh.ref_id()) |
160 |
+ |
{ |
161 |
+ |
if (lh.antisense_align() == rh.antisense_align()) |
162 |
+ |
fusion = true; |
163 |
+ |
else |
164 |
+ |
{ |
165 |
+ |
int inter_dist = 0; |
166 |
+ |
if (lh.antisense_align()) |
167 |
+ |
inter_dist = lh.left() - rh.right(); |
168 |
+ |
else |
169 |
+ |
inter_dist = rh.left() - lh.right(); |
170 |
+ |
|
171 |
+ |
if (inter_dist < -(int)max_insertion_length || inter_dist > (int)fusion_min_dist) |
172 |
+ |
fusion = true; |
173 |
+ |
} |
174 |
+ |
} |
175 |
+ |
} |
176 |
+ |
|
177 |
+ |
grade = InsertAlignmentGrade(lh, rh, min_mate_inner_dist, max_mate_inner_dist, fusion); |
178 |
+ |
} |
179 |
+ |
|
180 |
+ |
struct cmp_pair_alignment |
181 |
+ |
{ |
182 |
+ |
cmp_pair_alignment(const JunctionSet& gtf_juncs) : |
183 |
+ |
_gtf_juncs(>f_juncs) {} |
184 |
+ |
|
185 |
+ |
bool operator() (const pair<BowtieHit, BowtieHit>& l, const pair<BowtieHit, BowtieHit>& r) const |
186 |
+ |
{ |
187 |
+ |
InsertAlignmentGrade gl; set_insert_alignment_grade(l.first, l.second, gl); |
188 |
+ |
InsertAlignmentGrade gr; set_insert_alignment_grade(r.first, r.second, gr); |
189 |
+ |
|
190 |
+ |
bool better = gr < gl; |
191 |
+ |
bool worse = gl < gr; |
192 |
+ |
|
193 |
+ |
if (better && !worse) |
194 |
+ |
return true; |
195 |
+ |
else |
196 |
+ |
return false; |
197 |
+ |
} |
198 |
+ |
|
199 |
+ |
const JunctionSet* _gtf_juncs; |
200 |
+ |
}; |
201 |
+ |
|
202 |
+ |
void pair_best_alignments(const HitsForRead& left_hits, |
203 |
+ |
const HitsForRead& right_hits, |
204 |
+ |
InsertAlignmentGrade& best_grade, |
205 |
+ |
vector<pair<BowtieHit, BowtieHit> >& best_hits, |
206 |
+ |
const JunctionSet& gtf_junctions) |
207 |
+ |
{ |
208 |
|
const vector<BowtieHit>& left = left_hits.hits; |
209 |
|
const vector<BowtieHit>& right = right_hits.hits; |
210 |
|
|
217 |
|
const BowtieHit& rh = right[j]; |
218 |
|
if (rh.edit_dist() > max_read_mismatches) continue; |
219 |
|
|
220 |
< |
bool fusion = false; |
127 |
< |
if (fusion_search) |
128 |
< |
{ |
129 |
< |
bool left_fusion = lh.fusion_opcode() != FUSION_NOTHING; |
130 |
< |
bool right_fusion = rh.fusion_opcode() != FUSION_NOTHING; |
131 |
< |
if (left_fusion && right_fusion) |
132 |
< |
continue; |
133 |
< |
|
134 |
< |
fusion = left_fusion || right_fusion; |
135 |
< |
if (!fusion && lh.ref_id() != rh.ref_id()) |
136 |
< |
fusion = true; |
137 |
< |
|
138 |
< |
if (!fusion && lh.ref_id() == rh.ref_id()) |
139 |
< |
{ |
140 |
< |
if (lh.antisense_align() == rh.antisense_align()) |
141 |
< |
fusion = true; |
142 |
< |
else |
143 |
< |
{ |
144 |
< |
int inter_dist = 0; |
145 |
< |
if (lh.antisense_align()) |
146 |
< |
inter_dist = lh.left() - rh.right(); |
147 |
< |
else |
148 |
< |
inter_dist = rh.left() - lh.right(); |
149 |
< |
|
150 |
< |
if (inter_dist < -(int)max_insertion_length || inter_dist > (int)fusion_min_dist) |
151 |
< |
fusion = true; |
152 |
< |
} |
153 |
< |
} |
154 |
< |
} |
155 |
< |
|
156 |
< |
InsertAlignmentGrade g(lh, rh, min_mate_inner_dist, max_mate_inner_dist, fusion); |
220 |
> |
InsertAlignmentGrade g; set_insert_alignment_grade(lh, rh, g); |
221 |
|
|
222 |
< |
// Is the new status better than the current best one? |
159 |
< |
if (best_grade < g) |
222 |
> |
if (report_secondary_alignments) |
223 |
|
{ |
224 |
< |
left_best_hits.hits.clear(); |
162 |
< |
right_best_hits.hits.clear(); |
163 |
< |
|
164 |
< |
best_grade = g; |
165 |
< |
left_best_hits.hits.push_back(lh); |
166 |
< |
right_best_hits.hits.push_back(rh); |
224 |
> |
best_hits.push_back(make_pair(lh, rh)); |
225 |
|
} |
226 |
< |
else if (!(g < best_grade)) |
226 |
> |
else |
227 |
|
{ |
228 |
< |
left_best_hits.hits.push_back(lh); |
229 |
< |
right_best_hits.hits.push_back(rh); |
230 |
< |
best_grade.num_alignments++; |
228 |
> |
// Is the new status better than the current best one? |
229 |
> |
if (best_grade < g) |
230 |
> |
{ |
231 |
> |
best_hits.clear(); |
232 |
> |
best_grade = g; |
233 |
> |
best_hits.push_back(make_pair(lh, rh)); |
234 |
> |
} |
235 |
> |
else if (!(g < best_grade)) |
236 |
> |
{ |
237 |
> |
best_hits.push_back(make_pair(lh, rh)); |
238 |
> |
best_grade.num_alignments++; |
239 |
> |
} |
240 |
|
} |
241 |
|
} |
242 |
|
} |
243 |
+ |
|
244 |
+ |
if (report_secondary_alignments && best_hits.size() > 0) |
245 |
+ |
{ |
246 |
+ |
cmp_pair_alignment cmp(gtf_junctions); |
247 |
+ |
sort(best_hits.begin(), best_hits.end(), cmp); |
248 |
+ |
set_insert_alignment_grade(best_hits[0].first, best_hits[0].second, best_grade); |
249 |
+ |
best_grade.num_alignments = best_hits.size(); |
250 |
+ |
} |
251 |
|
} |
252 |
|
|
253 |
|
enum FragmentType {FRAG_UNPAIRED, FRAG_LEFT, FRAG_RIGHT}; |
255 |
|
void add_auxData(vector<string>& auxdata, vector<string>& sam_toks, |
256 |
|
const RefSequenceTable& rt,const BowtieHit& bh, FragmentType insert_side, |
257 |
|
int num_hits, const BowtieHit* next_hit, int hitIndex) { |
258 |
+ |
bool XS_found = false; |
259 |
+ |
if (sam_toks.size()>11) { |
260 |
|
|
261 |
< |
if (sam_toks.size()>11) { |
262 |
< |
for (size_t i=11;i<sam_toks.size();++i) { |
263 |
< |
if (sam_toks[i].find("NH:i:")==string::npos) |
264 |
< |
auxdata.push_back(sam_toks[i]); |
265 |
< |
} |
266 |
< |
} |
267 |
< |
string aux("NH:i:"); |
268 |
< |
str_appendInt(aux, num_hits); |
261 |
> |
for (size_t i=11;i<sam_toks.size();++i) { |
262 |
> |
if (sam_toks[i].find("NH:i:")==string::npos && |
263 |
> |
sam_toks[i].find("XS:i:")==string::npos) |
264 |
> |
auxdata.push_back(sam_toks[i]); |
265 |
> |
|
266 |
> |
if (sam_toks[i].find("XS:A:")!=string::npos) |
267 |
> |
XS_found = true; |
268 |
> |
} |
269 |
> |
} |
270 |
> |
string aux("NH:i:"); |
271 |
> |
str_appendInt(aux, num_hits); |
272 |
> |
auxdata.push_back(aux); |
273 |
> |
if (next_hit) { |
274 |
> |
const char* nh_ref_name = "="; |
275 |
> |
nh_ref_name = rt.get_name(next_hit->ref_id()); |
276 |
> |
assert (nh_ref_name != NULL); |
277 |
> |
bool same_contig=(next_hit->ref_id()==bh.ref_id()); |
278 |
> |
aux="CC:Z:"; |
279 |
> |
aux+= (same_contig ? "=" : nh_ref_name); |
280 |
> |
auxdata.push_back(aux); |
281 |
> |
aux="CP:i:"; |
282 |
> |
int nh_gpos=next_hit->left() + 1; |
283 |
> |
str_appendInt(aux, nh_gpos); |
284 |
|
auxdata.push_back(aux); |
285 |
< |
if (next_hit) { |
194 |
< |
const char* nh_ref_name = "="; |
195 |
< |
nh_ref_name = rt.get_name(next_hit->ref_id()); |
196 |
< |
assert (nh_ref_name != NULL); |
197 |
< |
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 |
285 |
> |
} //has next_hit |
286 |
|
// FIXME: this code is still a bit brittle, because it contains no |
287 |
|
// consistency check that the mates are on opposite strands - a current protocol |
288 |
|
// requirement, and that the strand indicated by the alignment is consistent |
289 |
|
// with the orientation of the splices (though that should be handled upstream). |
290 |
+ |
if (!XS_found) { |
291 |
|
const string xs_f("XS:A:+"); |
292 |
|
const string xs_r("XS:A:-"); |
293 |
< |
if (bh.contiguous()) { |
294 |
< |
if (library_type == FR_FIRSTSTRAND) { |
295 |
< |
if (insert_side == FRAG_LEFT || insert_side == FRAG_UNPAIRED) { |
296 |
< |
if (bh.antisense_align()) |
297 |
< |
auxdata.push_back(xs_f); |
298 |
< |
else |
299 |
< |
auxdata.push_back(xs_r); |
300 |
< |
} |
301 |
< |
else { |
302 |
< |
if (bh.antisense_align()) |
303 |
< |
auxdata.push_back(xs_r); |
304 |
< |
else |
305 |
< |
auxdata.push_back(xs_f); |
306 |
< |
} |
307 |
< |
} |
308 |
< |
else if (library_type == FR_SECONDSTRAND) { |
309 |
< |
if (insert_side == FRAG_LEFT || insert_side == FRAG_UNPAIRED){ |
310 |
< |
if (bh.antisense_align()) |
311 |
< |
auxdata.push_back(xs_r); |
312 |
< |
else |
313 |
< |
auxdata.push_back(xs_f); |
314 |
< |
} |
315 |
< |
else |
316 |
< |
{ |
317 |
< |
if (bh.antisense_align()) |
318 |
< |
auxdata.push_back(xs_f); |
319 |
< |
else |
320 |
< |
auxdata.push_back(xs_r); |
321 |
< |
} |
322 |
< |
} |
323 |
< |
} //bh.contiguous() |
243 |
< |
if (hitIndex >= 0) |
293 |
> |
if (library_type == FR_FIRSTSTRAND) { |
294 |
> |
if (insert_side == FRAG_LEFT || insert_side == FRAG_UNPAIRED) { |
295 |
> |
if (bh.antisense_align()) |
296 |
> |
auxdata.push_back(xs_f); |
297 |
> |
else |
298 |
> |
auxdata.push_back(xs_r); |
299 |
> |
} |
300 |
> |
else { |
301 |
> |
if (bh.antisense_align()) |
302 |
> |
auxdata.push_back(xs_r); |
303 |
> |
else |
304 |
> |
auxdata.push_back(xs_f); |
305 |
> |
} |
306 |
> |
} |
307 |
> |
else if (library_type == FR_SECONDSTRAND) { |
308 |
> |
if (insert_side == FRAG_LEFT || insert_side == FRAG_UNPAIRED){ |
309 |
> |
if (bh.antisense_align()) |
310 |
> |
auxdata.push_back(xs_r); |
311 |
> |
else |
312 |
> |
auxdata.push_back(xs_f); |
313 |
> |
} |
314 |
> |
else |
315 |
> |
{ |
316 |
> |
if (bh.antisense_align()) |
317 |
> |
auxdata.push_back(xs_f); |
318 |
> |
else |
319 |
> |
auxdata.push_back(xs_r); |
320 |
> |
} |
321 |
> |
} |
322 |
> |
} |
323 |
> |
if (hitIndex >= 0) |
324 |
|
{ |
325 |
< |
string aux("HI:i:"); |
326 |
< |
str_appendInt(aux, hitIndex); |
327 |
< |
auxdata.push_back(aux); |
325 |
> |
string aux("HI:i:"); |
326 |
> |
str_appendInt(aux, hitIndex); |
327 |
> |
auxdata.push_back(aux); |
328 |
|
} |
329 |
|
} |
330 |
|
|
599 |
|
|
600 |
|
void print_sam_for_single(const RefSequenceTable& rt, |
601 |
|
const HitsForRead& hits, |
602 |
< |
const FragmentAlignmentGrade& grade, |
603 |
< |
FragmentType frag_type, |
604 |
< |
//FILE* reads_file, |
605 |
< |
Read& read, |
606 |
< |
GBamWriter& bam_writer, |
607 |
< |
FILE* um_out//write unmapped reads here |
528 |
< |
) |
602 |
> |
const FragmentAlignmentGrade& grade, |
603 |
> |
FragmentType frag_type, |
604 |
> |
Read& read, |
605 |
> |
GBamWriter& bam_writer, |
606 |
> |
FILE* um_out //write unmapped reads here |
607 |
> |
) |
608 |
|
{ |
609 |
|
assert(!read.alt_name.empty()); |
610 |
+ |
if (hits.hits.empty()) |
611 |
+ |
return; |
612 |
+ |
|
613 |
|
lex_hit_sort s(rt, hits); |
614 |
|
vector<uint32_t> index_vector; |
615 |
|
|
619 |
|
|
620 |
|
sort(index_vector.begin(), index_vector.end(), s); |
621 |
|
|
622 |
< |
size_t primaryHit = (hits.hits.empty()? 0: random() % hits.hits.size()); |
622 |
> |
size_t primaryHit = 0; |
623 |
> |
if (!report_secondary_alignments) |
624 |
> |
primaryHit = random() % hits.hits.size(); |
625 |
> |
|
626 |
|
bool multipleHits = (hits.hits.size() > 1); |
627 |
|
for (i = 0; i < hits.hits.size(); ++i) |
628 |
|
{ |
544 |
– |
bool primary = (i == primaryHit); |
629 |
|
size_t index = index_vector[i]; |
630 |
+ |
bool primary = (index == primaryHit); |
631 |
|
const BowtieHit& bh = hits.hits[index]; |
632 |
|
rewrite_sam_record(bam_writer, rt, |
633 |
|
bh, |
643 |
|
} |
644 |
|
|
645 |
|
void print_sam_for_pair(const RefSequenceTable& rt, |
646 |
< |
const HitsForRead& left_hits, |
562 |
< |
const HitsForRead& right_hits, |
646 |
> |
const vector<pair<BowtieHit, BowtieHit> >& best_hits, |
647 |
|
const InsertAlignmentGrade& grade, |
648 |
|
ReadStream& left_reads_file, |
649 |
|
ReadStream& right_reads_file, |
653 |
|
uint64_t begin_id = 0, |
654 |
|
uint64_t end_id = std::numeric_limits<uint64_t>::max()) |
655 |
|
{ |
572 |
– |
assert (left_hits.insert_id == right_hits.insert_id); |
573 |
– |
|
656 |
|
Read left_read; |
657 |
|
Read right_read; |
658 |
< |
assert (left_hits.hits.size() == right_hits.hits.size() || |
659 |
< |
(left_hits.hits.empty() || right_hits.hits.empty())); |
658 |
> |
if (best_hits.empty()) |
659 |
> |
return; |
660 |
> |
|
661 |
> |
size_t i; |
662 |
> |
HitsForRead right_hits; |
663 |
> |
for (i = 0; i < best_hits.size(); ++i) |
664 |
> |
right_hits.hits.push_back(best_hits[i].second); |
665 |
|
|
666 |
|
size_t primaryHit = 0; |
667 |
|
vector<uint32_t> index_vector; |
668 |
< |
if(right_hits.hits.size() > 0) |
669 |
< |
{ |
670 |
< |
lex_hit_sort s(rt, right_hits); |
584 |
< |
for (size_t i = 0; i < right_hits.hits.size(); ++i) |
585 |
< |
index_vector.push_back(i); |
586 |
< |
|
587 |
< |
sort(index_vector.begin(), index_vector.end(), s); |
588 |
< |
primaryHit = random() % right_hits.hits.size(); |
589 |
< |
} |
590 |
< |
else if (left_hits.hits.size() > 0) |
591 |
< |
{ |
592 |
< |
lex_hit_sort s(rt, left_hits); |
593 |
< |
for (size_t i = 0; i < left_hits.hits.size(); ++i) |
594 |
< |
index_vector.push_back(i); |
595 |
< |
|
596 |
< |
sort(index_vector.begin(), index_vector.end(), s); |
597 |
< |
primaryHit = random() % left_hits.hits.size(); |
598 |
< |
} |
668 |
> |
lex_hit_sort s(rt, right_hits); |
669 |
> |
for (i = 0; i < right_hits.hits.size(); ++i) |
670 |
> |
index_vector.push_back(i); |
671 |
|
|
672 |
< |
bool got_left_read = left_reads_file.getRead(left_hits.insert_id, left_read, |
601 |
< |
reads_format, false, left_um_out, |
602 |
< |
false, begin_id, end_id); |
603 |
< |
|
604 |
< |
bool got_right_read = right_reads_file.getRead(right_hits.insert_id, right_read, |
605 |
< |
reads_format, false, right_um_out, |
606 |
< |
false, begin_id, end_id); |
607 |
< |
|
608 |
< |
if (left_hits.hits.size() == right_hits.hits.size()) |
609 |
< |
{ |
610 |
< |
assert (got_left_read && got_right_read); |
611 |
< |
bool multipleHits = (left_hits.hits.size() > 1); |
612 |
< |
for (size_t i = 0; i < right_hits.hits.size(); ++i) |
613 |
< |
{ |
614 |
< |
bool primary = (i == primaryHit); |
615 |
< |
size_t index = index_vector[i]; |
616 |
< |
const BowtieHit& right_bh = right_hits.hits[index]; |
617 |
< |
const BowtieHit& left_bh = left_hits.hits[index]; |
618 |
< |
|
619 |
< |
rewrite_sam_record(bam_writer, rt, |
620 |
< |
right_bh, |
621 |
< |
right_bh.hitfile_rec().c_str(), |
622 |
< |
right_read.alt_name.c_str(), |
623 |
< |
grade, |
624 |
< |
FRAG_RIGHT, |
625 |
< |
&left_bh, |
626 |
< |
right_hits.hits.size(), |
627 |
< |
(i < right_hits.hits.size() - 1) ? &(right_hits.hits[index_vector[i+1]]) : NULL, |
628 |
< |
primary, |
629 |
< |
(multipleHits? i: -1)); |
630 |
< |
rewrite_sam_record(bam_writer, rt, |
631 |
< |
left_bh, |
632 |
< |
left_bh.hitfile_rec().c_str(), |
633 |
< |
left_read.alt_name.c_str(), |
634 |
< |
grade, |
635 |
< |
FRAG_LEFT, |
636 |
< |
&right_bh, |
637 |
< |
left_hits.hits.size(), |
638 |
< |
(i < left_hits.hits.size() - 1) ? &(left_hits.hits[index_vector[i+1]]) : NULL, |
639 |
< |
primary, |
640 |
< |
(multipleHits? i: -1)); |
641 |
< |
} |
642 |
< |
} |
643 |
< |
else if (left_hits.hits.empty()) |
644 |
< |
{ //only right read was mapped properly |
645 |
< |
if (right_um_out) { |
646 |
< |
//write it in the mapped file with the #MAPPED# flag |
647 |
< |
fprintf(right_um_out, "@%s #MAPPED#\n%s\n+\n%s\n", right_read.alt_name.c_str(), |
648 |
< |
right_read.seq.c_str(), right_read.qual.c_str()); |
649 |
< |
} |
650 |
< |
for (size_t i = 0; i < right_hits.hits.size(); ++i) |
651 |
< |
{ |
652 |
< |
bool primary = (i == primaryHit); |
653 |
< |
size_t index = index_vector[i]; |
654 |
< |
const BowtieHit& bh = right_hits.hits[index]; |
655 |
< |
|
656 |
< |
rewrite_sam_record(bam_writer, rt, |
657 |
< |
bh, |
658 |
< |
bh.hitfile_rec().c_str(), |
659 |
< |
right_read.alt_name.c_str(), |
660 |
< |
grade, |
661 |
< |
FRAG_RIGHT, |
662 |
< |
NULL, |
663 |
< |
right_hits.hits.size(), |
664 |
< |
(i < right_hits.hits.size() - 1) ? &(right_hits.hits[index_vector[i+1]]) : NULL, |
665 |
< |
primary, |
666 |
< |
-1); |
667 |
< |
|
668 |
< |
} |
669 |
< |
} |
670 |
< |
else if (right_hits.hits.empty()) |
671 |
< |
{ //only left read was mapped properly |
672 |
< |
if (left_um_out) { |
673 |
< |
//write it in the mapped file with the #MAPPED# flag |
674 |
< |
fprintf(left_um_out, "@%s #MAPPED#\n%s\n+\n%s\n", left_read.alt_name.c_str(), left_read.seq.c_str(), |
675 |
< |
left_read.qual.c_str()); |
676 |
< |
} |
672 |
> |
sort(index_vector.begin(), index_vector.end(), s); |
673 |
|
|
674 |
< |
for (size_t i = 0; i < left_hits.hits.size(); ++i) |
675 |
< |
{ |
676 |
< |
bool primary = (i == primaryHit); |
677 |
< |
size_t index = index_vector[i]; |
678 |
< |
const BowtieHit& bh = left_hits.hits[index]; |
679 |
< |
rewrite_sam_record(bam_writer, rt, |
680 |
< |
bh, |
681 |
< |
bh.hitfile_rec().c_str(), |
682 |
< |
left_read.alt_name.c_str(), |
683 |
< |
grade, |
684 |
< |
FRAG_LEFT, |
685 |
< |
NULL, |
686 |
< |
left_hits.hits.size(), |
687 |
< |
(i < left_hits.hits.size() - 1) ? &(left_hits.hits[index_vector[i+1]]) : NULL, |
688 |
< |
primary, |
689 |
< |
-1); |
690 |
< |
|
691 |
< |
} |
692 |
< |
} |
693 |
< |
else |
694 |
< |
{ |
695 |
< |
assert (false); |
696 |
< |
} |
674 |
> |
if (!report_secondary_alignments) |
675 |
> |
primaryHit = random() % right_hits.hits.size(); |
676 |
> |
|
677 |
> |
bool got_left_read = left_reads_file.getRead(best_hits[0].first.insert_id(), left_read, |
678 |
> |
reads_format, false, begin_id, end_id, |
679 |
> |
left_um_out, false); |
680 |
> |
|
681 |
> |
bool got_right_read = right_reads_file.getRead(best_hits[0].first.insert_id(), right_read, |
682 |
> |
reads_format, false, begin_id, end_id, |
683 |
> |
right_um_out, false); |
684 |
> |
|
685 |
> |
assert (got_left_read && got_right_read); |
686 |
> |
bool multipleHits = (best_hits.size() > 1); |
687 |
> |
for (i = 0; i < best_hits.size(); ++i) |
688 |
> |
{ |
689 |
> |
size_t index = index_vector[i]; |
690 |
> |
bool primary = (index == primaryHit); |
691 |
> |
const BowtieHit& right_bh = best_hits[index].second; |
692 |
> |
const BowtieHit& left_bh = best_hits[index].first; |
693 |
> |
|
694 |
> |
rewrite_sam_record(bam_writer, rt, |
695 |
> |
right_bh, |
696 |
> |
right_bh.hitfile_rec().c_str(), |
697 |
> |
right_read.alt_name.c_str(), |
698 |
> |
grade, |
699 |
> |
FRAG_RIGHT, |
700 |
> |
&left_bh, |
701 |
> |
best_hits.size(), |
702 |
> |
(i < best_hits.size() - 1) ? &(best_hits[index_vector[i+1]].second) : NULL, |
703 |
> |
primary, |
704 |
> |
(multipleHits? i: -1)); |
705 |
> |
rewrite_sam_record(bam_writer, rt, |
706 |
> |
left_bh, |
707 |
> |
left_bh.hitfile_rec().c_str(), |
708 |
> |
left_read.alt_name.c_str(), |
709 |
> |
grade, |
710 |
> |
FRAG_LEFT, |
711 |
> |
&right_bh, |
712 |
> |
best_hits.size(), |
713 |
> |
(i < best_hits.size() - 1) ? &(best_hits[index_vector[i+1]].first) : NULL, |
714 |
> |
primary, |
715 |
> |
(multipleHits? i: -1)); |
716 |
> |
} |
717 |
|
} |
718 |
|
|
719 |
|
/** |
801 |
|
for (size_t i = 0; i < hits.hits.size(); ++i) |
802 |
|
{ |
803 |
|
BowtieHit& bh = hits.hits[i]; |
804 |
+ |
if (bh.edit_dist() > max_read_mismatches) |
805 |
+ |
continue; |
806 |
+ |
|
807 |
|
bool filter_hit = false; |
808 |
|
if (!bh.contiguous()) |
809 |
|
{ |
865 |
|
{ |
866 |
|
HitsForRead best_hits; |
867 |
|
best_hits.insert_id = curr_left_obs_order; |
868 |
< |
FragmentAlignmentGrade grade; |
868 |
> |
FragmentAlignmentGrade best_grade; |
869 |
|
|
870 |
|
// Process hits for left singleton, select best alignments |
871 |
< |
read_best_alignments(curr_left_hit_group, grade, best_hits, *gtf_junctions); |
871 |
> |
read_best_alignments(curr_left_hit_group, best_grade, best_hits, *gtf_junctions); |
872 |
|
update_junctions(best_hits, *junctions); |
873 |
|
update_fusions(best_hits, *rt, *fusions); |
874 |
|
|
884 |
|
{ |
885 |
|
HitsForRead best_hits; |
886 |
|
best_hits.insert_id = curr_right_obs_order; |
887 |
+ |
FragmentAlignmentGrade best_grade; |
888 |
|
if (curr_right_obs_order >= begin_id) |
889 |
|
{ |
870 |
– |
FragmentAlignmentGrade grade; |
890 |
|
// Process hit for right singleton, select best alignments |
891 |
< |
read_best_alignments(curr_right_hit_group,grade, best_hits, *gtf_junctions); |
891 |
> |
read_best_alignments(curr_right_hit_group, best_grade, best_hits, *gtf_junctions); |
892 |
|
update_junctions(best_hits, *junctions); |
893 |
|
update_fusions(best_hits, *rt, *fusions); |
894 |
|
} |
904 |
|
curr_left_obs_order < end_id && |
905 |
|
curr_left_obs_order != VMAXINT32) |
906 |
|
{ |
907 |
< |
if (curr_left_hit_group.hits.empty()) |
908 |
< |
{ |
909 |
< |
HitsForRead right_best_hits; |
910 |
< |
right_best_hits.insert_id = curr_right_obs_order; |
911 |
< |
|
912 |
< |
FragmentAlignmentGrade grade; |
913 |
< |
read_best_alignments(curr_right_hit_group, grade, right_best_hits, *gtf_junctions); |
914 |
< |
update_junctions(right_best_hits, *junctions); |
915 |
< |
update_fusions(right_best_hits, *rt, *fusions); |
916 |
< |
} |
917 |
< |
else if (curr_right_hit_group.hits.empty()) |
907 |
> |
vector<pair<BowtieHit, BowtieHit> > best_hits; |
908 |
> |
|
909 |
> |
InsertAlignmentGrade grade; |
910 |
> |
pair_best_alignments(curr_left_hit_group, |
911 |
> |
curr_right_hit_group, |
912 |
> |
grade, |
913 |
> |
best_hits, |
914 |
> |
*gtf_junctions); |
915 |
> |
|
916 |
> |
HitsForRead single_best_hits; |
917 |
> |
single_best_hits.insert_id = curr_left_obs_order; |
918 |
> |
|
919 |
> |
// this is too much memory copy - we will make it efficient soon |
920 |
> |
for (size_t i = 0; i < best_hits.size(); ++i) |
921 |
|
{ |
922 |
< |
HitsForRead left_best_hits; |
923 |
< |
left_best_hits.insert_id = curr_left_obs_order; |
902 |
< |
|
903 |
< |
FragmentAlignmentGrade grade; |
904 |
< |
// Process hits for left singleton, select best alignments |
905 |
< |
read_best_alignments(curr_left_hit_group, grade, left_best_hits, *gtf_junctions); |
906 |
< |
update_junctions(left_best_hits, *junctions); |
907 |
< |
update_fusions(left_best_hits, *rt, *fusions); |
922 |
> |
single_best_hits.hits.push_back(best_hits[i].first); |
923 |
> |
single_best_hits.hits.push_back(best_hits[i].second); |
924 |
|
} |
909 |
– |
else |
910 |
– |
{ |
911 |
– |
HitsForRead left_best_hits; |
912 |
– |
HitsForRead right_best_hits; |
913 |
– |
left_best_hits.insert_id = curr_left_obs_order; |
914 |
– |
right_best_hits.insert_id = curr_right_obs_order; |
915 |
– |
|
916 |
– |
// daehwan - apply gtf_junctions here, too! |
925 |
|
|
926 |
< |
InsertAlignmentGrade grade; |
927 |
< |
pair_best_alignments(curr_left_hit_group, |
928 |
< |
curr_right_hit_group, |
921 |
< |
grade, |
922 |
< |
left_best_hits, |
923 |
< |
right_best_hits); |
924 |
< |
update_junctions(left_best_hits, *junctions); |
925 |
< |
update_junctions(right_best_hits, *junctions); |
926 |
< |
update_fusions(left_best_hits, *rt, *fusions); |
927 |
< |
update_fusions(right_best_hits, *rt, *fusions); |
928 |
< |
} |
929 |
< |
|
926 |
> |
update_junctions(single_best_hits, *junctions); |
927 |
> |
update_fusions(single_best_hits, *rt, *fusions); |
928 |
> |
|
929 |
|
l_hs.next_read_hits(curr_left_hit_group); |
930 |
|
curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id); |
931 |
|
|
1020 |
|
best_hits.insert_id = curr_left_obs_order; |
1021 |
|
FragmentAlignmentGrade grade; |
1022 |
|
bool got_read=left_reads_file.getRead(curr_left_obs_order, l_read, |
1023 |
< |
reads_format, false, left_um_out.file, |
1024 |
< |
false, begin_id, end_id); |
1025 |
< |
assert(got_read); |
1023 |
> |
reads_format, false, begin_id, end_id, |
1024 |
> |
left_um_out.file, false); |
1025 |
> |
//assert(got_read); |
1026 |
|
|
1027 |
|
if (right_reads_file.file()) { |
1028 |
+ |
/* |
1029 |
|
fprintf(left_um_out.file, "@%s #MAPPED#\n%s\n+\n%s\n", l_read.alt_name.c_str(), |
1030 |
|
l_read.seq.c_str(), l_read.qual.c_str()); |
1031 |
+ |
*/ |
1032 |
|
got_read=right_reads_file.getRead(curr_left_obs_order, r_read, |
1033 |
< |
reads_format, false, right_um_out.file, |
1034 |
< |
true, begin_id, end_id); |
1035 |
< |
assert(got_read); |
1033 |
> |
reads_format, false, begin_id, end_id, |
1034 |
> |
right_um_out.file, true); |
1035 |
> |
//assert(got_read); |
1036 |
|
} |
1037 |
< |
|
1037 |
> |
|
1038 |
|
exclude_hits_on_filtered_junctions(*junctions, curr_left_hit_group); |
1039 |
|
|
1040 |
|
// Process hits for left singleton, select best alignments |
1041 |
|
read_best_alignments(curr_left_hit_group, grade, best_hits, *gtf_junctions); |
1042 |
|
|
1043 |
< |
if (best_hits.hits.size()>0 && best_hits.hits.size() <= max_multihits) |
1043 |
> |
if (best_hits.hits.size() > max_multihits) |
1044 |
> |
{ |
1045 |
> |
best_hits.hits.erase(best_hits.hits.begin() + max_multihits, best_hits.hits.end()); |
1046 |
> |
grade.num_alignments = best_hits.hits.size(); |
1047 |
> |
} |
1048 |
> |
|
1049 |
> |
if (best_hits.hits.size() > 0) |
1050 |
|
{ |
1051 |
|
update_junctions(best_hits, *final_junctions); |
1052 |
|
update_insertions_and_deletions(best_hits, *final_insertions, *final_deletions); |
1077 |
|
if (curr_right_obs_order >= begin_id) |
1078 |
|
{ |
1079 |
|
bool got_read=right_reads_file.getRead(curr_right_obs_order, r_read, |
1080 |
< |
reads_format, false, right_um_out.file, |
1081 |
< |
false, begin_id, end_id); |
1075 |
< |
|
1076 |
< |
assert(got_read); |
1080 |
> |
reads_format, false, begin_id, end_id, |
1081 |
> |
right_um_out.file, false); |
1082 |
|
|
1083 |
+ |
//assert(got_read); |
1084 |
+ |
/* |
1085 |
|
fprintf(right_um_out.file, "@%s #MAPPED#\n%s\n+\n%s\n", r_read.alt_name.c_str(), |
1086 |
|
r_read.seq.c_str(), r_read.qual.c_str()); |
1087 |
+ |
*/ |
1088 |
|
got_read=left_reads_file.getRead(curr_right_obs_order, l_read, |
1089 |
< |
reads_format, false, left_um_out.file, |
1090 |
< |
true, begin_id, end_id); |
1091 |
< |
assert(got_read); |
1089 |
> |
reads_format, false, begin_id, end_id, |
1090 |
> |
left_um_out.file, true); |
1091 |
> |
//assert(got_read); |
1092 |
|
|
1093 |
|
exclude_hits_on_filtered_junctions(*junctions, curr_right_hit_group); |
1094 |
|
|
1095 |
|
// Process hit for right singleton, select best alignments |
1096 |
|
read_best_alignments(curr_right_hit_group, grade, best_hits, *gtf_junctions); |
1097 |
< |
if (best_hits.hits.size()>0 && best_hits.hits.size() <= max_multihits) |
1097 |
> |
|
1098 |
> |
if (best_hits.hits.size() > max_multihits) |
1099 |
> |
{ |
1100 |
> |
best_hits.hits.erase(best_hits.hits.begin() + max_multihits, best_hits.hits.end()); |
1101 |
> |
grade.num_alignments = best_hits.hits.size(); |
1102 |
> |
} |
1103 |
> |
|
1104 |
> |
if (best_hits.hits.size() > 0) |
1105 |
|
{ |
1106 |
|
update_junctions(best_hits, *final_junctions); |
1107 |
|
update_insertions_and_deletions(best_hits, *final_insertions, *final_deletions); |
1135 |
|
//write it in the mapped file with the #MAPPED# flag |
1136 |
|
|
1137 |
|
bool got_read=right_reads_file.getRead(curr_left_obs_order, r_read, |
1138 |
< |
reads_format, false, right_um_out.file, |
1139 |
< |
false, begin_id, end_id); |
1140 |
< |
assert(got_read); |
1138 |
> |
reads_format, false, begin_id, end_id, |
1139 |
> |
right_um_out.file, false); |
1140 |
> |
//assert(got_read); |
1141 |
> |
/* |
1142 |
|
fprintf(right_um_out.file, "@%s #MAPPED#\n%s\n+\n%s\n", r_read.alt_name.c_str(), |
1143 |
|
r_read.seq.c_str(), r_read.qual.c_str()); |
1144 |
+ |
*/ |
1145 |
|
HitsForRead right_best_hits; |
1146 |
|
right_best_hits.insert_id = curr_right_obs_order; |
1147 |
|
|
1148 |
|
FragmentAlignmentGrade grade; |
1149 |
|
read_best_alignments(curr_right_hit_group, grade, right_best_hits, *gtf_junctions); |
1150 |
+ |
|
1151 |
+ |
if (right_best_hits.hits.size() > max_multihits) |
1152 |
+ |
{ |
1153 |
+ |
right_best_hits.hits.erase(right_best_hits.hits.begin() + max_multihits, right_best_hits.hits.end()); |
1154 |
+ |
grade.num_alignments = right_best_hits.hits.size(); |
1155 |
+ |
} |
1156 |
|
|
1157 |
< |
if (right_best_hits.hits.size()>0 && right_best_hits.hits.size() <= max_multihits) |
1157 |
> |
if (right_best_hits.hits.size() > 0) |
1158 |
|
{ |
1159 |
|
update_junctions(right_best_hits, *final_junctions); |
1160 |
|
update_insertions_and_deletions(right_best_hits, *final_insertions, *final_deletions); |
1175 |
|
left_best_hits.insert_id = curr_left_obs_order; |
1176 |
|
//only left read mapped |
1177 |
|
bool got_read=left_reads_file.getRead(curr_left_obs_order, l_read, |
1178 |
< |
reads_format, false, left_um_out.file, |
1179 |
< |
false, begin_id, end_id); |
1180 |
< |
assert(got_read); |
1178 |
> |
reads_format, false, begin_id, end_id, |
1179 |
> |
left_um_out.file, false); |
1180 |
> |
//assert(got_read); |
1181 |
> |
/* |
1182 |
|
fprintf(left_um_out.file, "@%s #MAPPED#\n%s\n+\n%s\n", l_read.alt_name.c_str(), |
1183 |
|
l_read.seq.c_str(), l_read.qual.c_str()); |
1184 |
+ |
*/ |
1185 |
|
FragmentAlignmentGrade grade; |
1186 |
|
// Process hits for left singleton, select best alignments |
1187 |
|
read_best_alignments(curr_left_hit_group, grade, left_best_hits, *gtf_junctions); |
1188 |
< |
if (left_best_hits.hits.size()>0 && left_best_hits.hits.size() <= max_multihits) |
1188 |
> |
|
1189 |
> |
if (left_best_hits.hits.size() > max_multihits) |
1190 |
> |
{ |
1191 |
> |
left_best_hits.hits.erase(left_best_hits.hits.begin() + max_multihits, left_best_hits.hits.end()); |
1192 |
> |
grade.num_alignments = left_best_hits.hits.size(); |
1193 |
> |
} |
1194 |
> |
|
1195 |
> |
if (left_best_hits.hits.size() > 0) |
1196 |
|
{ |
1197 |
|
update_junctions(left_best_hits, *final_junctions); |
1198 |
|
update_insertions_and_deletions(left_best_hits, *final_insertions, *final_deletions); |
1208 |
|
} |
1209 |
|
} |
1210 |
|
else |
1211 |
< |
{ //hits for both left and right reads |
1212 |
< |
HitsForRead left_best_hits; |
1181 |
< |
HitsForRead right_best_hits; |
1182 |
< |
left_best_hits.insert_id = curr_left_obs_order; |
1183 |
< |
right_best_hits.insert_id = curr_right_obs_order; |
1211 |
> |
{ |
1212 |
> |
vector<pair<BowtieHit, BowtieHit> > best_hits; |
1213 |
|
|
1214 |
|
InsertAlignmentGrade grade; |
1215 |
|
pair_best_alignments(curr_left_hit_group, |
1216 |
|
curr_right_hit_group, |
1217 |
|
grade, |
1218 |
< |
left_best_hits, |
1219 |
< |
right_best_hits); |
1218 |
> |
best_hits, |
1219 |
> |
*gtf_junctions); |
1220 |
> |
|
1221 |
> |
if (best_hits.size() > max_multihits) |
1222 |
> |
{ |
1223 |
> |
best_hits.erase(best_hits.begin() + max_multihits, best_hits.end()); |
1224 |
> |
grade.num_alignments = best_hits.size(); |
1225 |
> |
} |
1226 |
> |
|
1227 |
> |
//hits for both left and right reads |
1228 |
> |
HitsForRead single_best_hits; |
1229 |
> |
single_best_hits.insert_id = curr_left_obs_order; |
1230 |
> |
for (size_t i = 0; i < best_hits.size(); ++i) |
1231 |
> |
{ |
1232 |
> |
single_best_hits.hits.push_back(best_hits[i].first); |
1233 |
> |
single_best_hits.hits.push_back(best_hits[i].second); |
1234 |
> |
} |
1235 |
|
|
1236 |
< |
if (left_best_hits.hits.size()>0 && left_best_hits.hits.size() <= max_multihits && |
1193 |
< |
right_best_hits.hits.size()>0 && right_best_hits.hits.size() <= max_multihits) |
1236 |
> |
if (single_best_hits.hits.size() > 0) |
1237 |
|
{ |
1238 |
< |
update_junctions(left_best_hits, *final_junctions); |
1239 |
< |
update_junctions(right_best_hits, *final_junctions); |
1197 |
< |
update_insertions_and_deletions(left_best_hits, *final_insertions, *final_deletions); |
1198 |
< |
update_insertions_and_deletions(right_best_hits, *final_insertions, *final_deletions); |
1238 |
> |
update_junctions(single_best_hits, *final_junctions); |
1239 |
> |
update_insertions_and_deletions(single_best_hits, *final_insertions, *final_deletions); |
1240 |
|
|
1241 |
< |
pair_support(left_best_hits, right_best_hits, *final_fusions, *fusions); |
1242 |
< |
update_fusions(left_best_hits, *rt, *final_fusions, *fusions); |
1202 |
< |
update_fusions(right_best_hits, *rt, *final_fusions, *fusions); |
1241 |
> |
pair_support(best_hits, *final_fusions, *fusions); |
1242 |
> |
update_fusions(single_best_hits, *rt, *final_fusions, *fusions); |
1243 |
|
|
1244 |
|
print_sam_for_pair(*rt, |
1245 |
< |
left_best_hits, |
1206 |
< |
right_best_hits, |
1245 |
> |
best_hits, |
1246 |
|
grade, |
1247 |
|
left_reads_file, |
1248 |
|
right_reads_file, |
1263 |
|
} //while we still have unreported hits.. |
1264 |
|
//print the remaining unmapped reads at the end of each reads' stream |
1265 |
|
left_reads_file.getRead(VMAXINT32, l_read, |
1266 |
< |
reads_format, false, left_um_out.file, |
1267 |
< |
false, begin_id, end_id); |
1266 |
> |
reads_format, false, begin_id, end_id, |
1267 |
> |
left_um_out.file, false); |
1268 |
|
if (right_reads_file.file()) |
1269 |
|
right_reads_file.getRead(VMAXINT32, r_read, |
1270 |
< |
reads_format, false, right_um_out.file, |
1271 |
< |
false, begin_id, end_id); |
1270 |
> |
reads_format, false, begin_id, end_id, |
1271 |
> |
right_um_out.file, false); |
1272 |
|
|
1273 |
|
left_um_out.close(); |
1274 |
|
right_um_out.close(); |