44 |
|
#include "wiggles.h" |
45 |
|
#include "tokenize.h" |
46 |
|
#include "reads.h" |
47 |
+ |
#include "coverage.h" |
48 |
|
|
49 |
|
|
50 |
|
#include "inserts.h" |
53 |
|
using namespace seqan; |
54 |
|
using std::set; |
55 |
|
|
56 |
+ |
static const JunctionSet empty_junctions; |
57 |
+ |
static const InsertionSet empty_insertions; |
58 |
+ |
static const DeletionSet empty_deletions; |
59 |
+ |
static const FusionSet empty_fusions; |
60 |
+ |
static const Coverage empty_coverage; |
61 |
+ |
|
62 |
|
// daehwan - this is redundancy, which should be removed. |
63 |
|
void get_seqs(istream& ref_stream, |
64 |
|
RefSequenceTable& rt, |
80 |
|
} |
81 |
|
} |
82 |
|
|
83 |
+ |
struct cmp_read_alignment |
84 |
+ |
{ |
85 |
+ |
bool operator() (const BowtieHit& l, const BowtieHit& r) const |
86 |
+ |
{ |
87 |
+ |
return l.alignment_score() > r.alignment_score(); |
88 |
+ |
} |
89 |
+ |
}; |
90 |
+ |
|
91 |
|
void read_best_alignments(const HitsForRead& hits_for_read, |
92 |
< |
FragmentAlignmentGrade& best_grade, |
93 |
< |
HitsForRead& best_hits, |
94 |
< |
const JunctionSet& gtf_junctions) |
92 |
> |
HitsForRead& best_hits, |
93 |
> |
const JunctionSet& gtf_junctions, |
94 |
> |
const JunctionSet& junctions = empty_junctions, |
95 |
> |
const InsertionSet& insertions = empty_insertions, |
96 |
> |
const DeletionSet& deletions = empty_deletions, |
97 |
> |
const FusionSet& fusions = empty_fusions, |
98 |
> |
const Coverage& coverage = empty_coverage, |
99 |
> |
bool final_report = false) |
100 |
|
{ |
101 |
|
const vector<BowtieHit>& hits = hits_for_read.hits; |
102 |
+ |
|
103 |
+ |
if (hits.size() >= max_multihits * 5) |
104 |
+ |
return; |
105 |
+ |
|
106 |
|
for (size_t i = 0; i < hits.size(); ++i) |
107 |
|
{ |
108 |
< |
if (hits[i].edit_dist()>max_read_mismatches) continue; |
109 |
< |
|
110 |
< |
FragmentAlignmentGrade g(hits[i], gtf_junctions); |
111 |
< |
|
112 |
< |
// Is the new status better than the current best one? |
113 |
< |
if (best_grade < g) |
114 |
< |
{ |
115 |
< |
best_hits.hits.clear(); |
116 |
< |
best_grade = g; |
117 |
< |
best_hits.hits.push_back(hits[i]); |
118 |
< |
} |
119 |
< |
else if (!(g < best_grade)) // is it just as good? |
120 |
< |
{ |
121 |
< |
best_grade.num_alignments++; |
122 |
< |
best_hits.hits.push_back(hits[i]); |
123 |
< |
} |
108 |
> |
if (hits[i].edit_dist() > max_read_mismatches) |
109 |
> |
continue; |
110 |
> |
|
111 |
> |
BowtieHit hit = hits[i]; |
112 |
> |
AlignStatus align_status(hit, gtf_junctions, |
113 |
> |
junctions, insertions, deletions, fusions, coverage); |
114 |
> |
hit.alignment_score(align_status._alignment_score); |
115 |
> |
|
116 |
> |
if (report_secondary_alignments || !final_report) |
117 |
> |
{ |
118 |
> |
best_hits.hits.push_back(hit); |
119 |
> |
} |
120 |
> |
else |
121 |
> |
{ |
122 |
> |
// Is the new status better than the current best one? |
123 |
> |
if (best_hits.hits.size() == 0 || cmp_read_alignment()(hit, best_hits.hits[0])) |
124 |
> |
{ |
125 |
> |
best_hits.hits.clear(); |
126 |
> |
best_hits.hits.push_back(hit); |
127 |
> |
} |
128 |
> |
else if (!cmp_read_alignment()(best_hits.hits[0], hit)) // is it just as good? |
129 |
> |
{ |
130 |
> |
best_hits.hits.push_back(hit); |
131 |
> |
} |
132 |
> |
} |
133 |
> |
} |
134 |
> |
|
135 |
> |
if ((report_secondary_alignments || !final_report) && best_hits.hits.size() > 0) |
136 |
> |
{ |
137 |
> |
sort(best_hits.hits.begin(), best_hits.hits.end(), cmp_read_alignment()); |
138 |
> |
} |
139 |
> |
|
140 |
> |
if (final_report) |
141 |
> |
{ |
142 |
> |
if (best_hits.hits.size() > max_multihits) |
143 |
> |
best_hits.hits.erase(best_hits.hits.begin() + max_multihits, best_hits.hits.end()); |
144 |
|
} |
145 |
|
} |
146 |
|
|
147 |
< |
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) |
147 |
> |
bool set_insert_alignment_grade(const BowtieHit& lh, const BowtieHit& rh, InsertAlignmentGrade& grade) |
148 |
|
{ |
149 |
|
// max mate inner distance (genomic) |
150 |
|
int min_mate_inner_dist = inner_dist_mean - inner_dist_std_dev; |
151 |
|
if (max_mate_inner_dist == -1) |
152 |
|
max_mate_inner_dist = inner_dist_mean + inner_dist_std_dev; |
153 |
|
|
154 |
+ |
bool fusion = false; |
155 |
+ |
bool left_fusion = lh.fusion_opcode() != FUSION_NOTHING; |
156 |
+ |
bool right_fusion = rh.fusion_opcode() != FUSION_NOTHING; |
157 |
+ |
if (left_fusion && right_fusion) |
158 |
+ |
return false; |
159 |
+ |
|
160 |
+ |
fusion = left_fusion || right_fusion; |
161 |
+ |
if (!fusion && lh.ref_id() != rh.ref_id()) |
162 |
+ |
fusion = true; |
163 |
+ |
|
164 |
+ |
if (!fusion && lh.ref_id() == rh.ref_id()) |
165 |
+ |
{ |
166 |
+ |
if (lh.antisense_align() == rh.antisense_align()) |
167 |
+ |
fusion = true; |
168 |
+ |
else |
169 |
+ |
{ |
170 |
+ |
int inner_dist = 0; |
171 |
+ |
if (lh.antisense_align()) |
172 |
+ |
// used rh.left() instead of rh.right() for the cases, |
173 |
+ |
// where reads overlap with each other or reads span introns |
174 |
+ |
inner_dist = lh.left() - rh.left(); |
175 |
+ |
else |
176 |
+ |
inner_dist = rh.left() - lh.left(); |
177 |
+ |
|
178 |
+ |
if (inner_dist < 0 || inner_dist > (int)fusion_min_dist) |
179 |
+ |
fusion = true; |
180 |
+ |
} |
181 |
+ |
} |
182 |
+ |
|
183 |
+ |
grade = InsertAlignmentGrade(lh, rh, min_mate_inner_dist, max_mate_inner_dist, fusion); |
184 |
+ |
|
185 |
+ |
return true; |
186 |
+ |
} |
187 |
+ |
|
188 |
+ |
struct cmp_pair_alignment |
189 |
+ |
{ |
190 |
+ |
cmp_pair_alignment(const JunctionSet& gtf_juncs) : |
191 |
+ |
_gtf_juncs(>f_juncs) {} |
192 |
+ |
|
193 |
+ |
bool operator() (const pair<BowtieHit, BowtieHit>& l, const pair<BowtieHit, BowtieHit>& r) const |
194 |
+ |
{ |
195 |
+ |
InsertAlignmentGrade gl; set_insert_alignment_grade(l.first, l.second, gl); |
196 |
+ |
InsertAlignmentGrade gr; set_insert_alignment_grade(r.first, r.second, gr); |
197 |
+ |
|
198 |
+ |
bool better = gr < gl; |
199 |
+ |
bool worse = gl < gr; |
200 |
+ |
|
201 |
+ |
if (better && !worse) |
202 |
+ |
return true; |
203 |
+ |
else |
204 |
+ |
return false; |
205 |
+ |
} |
206 |
+ |
|
207 |
+ |
const JunctionSet* _gtf_juncs; |
208 |
+ |
}; |
209 |
+ |
|
210 |
+ |
void pair_best_alignments(const HitsForRead& left_hits, |
211 |
+ |
const HitsForRead& right_hits, |
212 |
+ |
InsertAlignmentGrade& best_grade, |
213 |
+ |
vector<pair<BowtieHit, BowtieHit> >& best_hits, |
214 |
+ |
const JunctionSet& gtf_junctions, |
215 |
+ |
const JunctionSet& junctions = empty_junctions, |
216 |
+ |
const InsertionSet& insertions = empty_insertions, |
217 |
+ |
const DeletionSet& deletions = empty_deletions, |
218 |
+ |
const FusionSet& fusions = empty_fusions, |
219 |
+ |
const Coverage& coverage = empty_coverage, |
220 |
+ |
bool final_report = false) |
221 |
+ |
{ |
222 |
|
const vector<BowtieHit>& left = left_hits.hits; |
223 |
|
const vector<BowtieHit>& right = right_hits.hits; |
224 |
< |
|
224 |
> |
|
225 |
> |
if (left.size() >= max_multihits * 5 || right.size() >= max_multihits * 5) |
226 |
> |
return; |
227 |
> |
|
228 |
|
for (size_t i = 0; i < left.size(); ++i) |
229 |
|
{ |
230 |
< |
const BowtieHit& lh = left[i]; |
231 |
< |
if (lh.edit_dist() > max_read_mismatches) continue; |
230 |
> |
if (left[i].edit_dist() > max_read_mismatches) continue; |
231 |
> |
|
232 |
> |
BowtieHit lh = left[i]; |
233 |
> |
AlignStatus align_status(lh, gtf_junctions, |
234 |
> |
junctions, insertions, deletions, fusions, coverage); |
235 |
> |
lh.alignment_score(align_status._alignment_score); |
236 |
> |
|
237 |
|
for (size_t j = 0; j < right.size(); ++j) |
238 |
|
{ |
239 |
< |
const BowtieHit& rh = right[j]; |
240 |
< |
if (rh.edit_dist() > max_read_mismatches) continue; |
241 |
< |
|
242 |
< |
bool fusion = false; |
243 |
< |
if (fusion_search) |
244 |
< |
{ |
245 |
< |
bool left_fusion = lh.fusion_opcode() != FUSION_NOTHING; |
246 |
< |
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 |
< |
} |
239 |
> |
if (right[j].edit_dist() > max_read_mismatches) continue; |
240 |
> |
|
241 |
> |
BowtieHit rh = right[j]; |
242 |
> |
AlignStatus align_status(rh, gtf_junctions, |
243 |
> |
junctions, insertions, deletions, fusions, coverage); |
244 |
> |
rh.alignment_score(align_status._alignment_score); |
245 |
> |
InsertAlignmentGrade g; |
246 |
> |
bool allowed = set_insert_alignment_grade(lh, rh, g); |
247 |
|
|
248 |
< |
InsertAlignmentGrade g(lh, rh, min_mate_inner_dist, max_mate_inner_dist, fusion); |
248 |
> |
if (!allowed) continue; |
249 |
> |
if (!fusion_search && !report_discordant_pair_alignments && g.is_fusion()) continue; |
250 |
|
|
251 |
< |
// Is the new status better than the current best one? |
159 |
< |
if (best_grade < g) |
251 |
> |
if (report_secondary_alignments || !final_report) |
252 |
|
{ |
253 |
< |
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); |
253 |
> |
best_hits.push_back(make_pair(lh, rh)); |
254 |
|
} |
255 |
< |
else if (!(g < best_grade)) |
255 |
> |
else |
256 |
|
{ |
257 |
< |
left_best_hits.hits.push_back(lh); |
258 |
< |
right_best_hits.hits.push_back(rh); |
259 |
< |
best_grade.num_alignments++; |
257 |
> |
// Is the new status better than the current best one? |
258 |
> |
if (best_grade < g) |
259 |
> |
{ |
260 |
> |
best_hits.clear(); |
261 |
> |
best_grade = g; |
262 |
> |
best_hits.push_back(make_pair(lh, rh)); |
263 |
> |
} |
264 |
> |
else if (!(g < best_grade)) |
265 |
> |
{ |
266 |
> |
best_hits.push_back(make_pair(lh, rh)); |
267 |
> |
} |
268 |
|
} |
269 |
|
} |
270 |
|
} |
271 |
+ |
|
272 |
+ |
if ((report_secondary_alignments || !final_report) && best_hits.size() > 0) |
273 |
+ |
{ |
274 |
+ |
cmp_pair_alignment cmp(gtf_junctions); |
275 |
+ |
sort(best_hits.begin(), best_hits.end(), cmp); |
276 |
+ |
set_insert_alignment_grade(best_hits[0].first, best_hits[0].second, best_grade); |
277 |
+ |
} |
278 |
+ |
|
279 |
+ |
if (final_report) |
280 |
+ |
{ |
281 |
+ |
if (best_hits.size() > max_multihits) |
282 |
+ |
best_hits.erase(best_hits.begin() + max_multihits, best_hits.end()); |
283 |
+ |
} |
284 |
+ |
|
285 |
+ |
best_grade.num_alignments = best_hits.size(); |
286 |
|
} |
287 |
|
|
288 |
|
enum FragmentType {FRAG_UNPAIRED, FRAG_LEFT, FRAG_RIGHT}; |
290 |
|
void add_auxData(vector<string>& auxdata, vector<string>& sam_toks, |
291 |
|
const RefSequenceTable& rt,const BowtieHit& bh, FragmentType insert_side, |
292 |
|
int num_hits, const BowtieHit* next_hit, int hitIndex) { |
293 |
+ |
bool XS_found = false; |
294 |
+ |
if (sam_toks.size()>11) { |
295 |
|
|
296 |
< |
if (sam_toks.size()>11) { |
297 |
< |
for (size_t i=11;i<sam_toks.size();++i) { |
298 |
< |
if (sam_toks[i].find("NH:i:")==string::npos) |
299 |
< |
auxdata.push_back(sam_toks[i]); |
300 |
< |
} |
301 |
< |
} |
302 |
< |
string aux("NH:i:"); |
303 |
< |
str_appendInt(aux, num_hits); |
296 |
> |
for (size_t i=11;i<sam_toks.size();++i) { |
297 |
> |
if (sam_toks[i].find("NH:i:")==string::npos && |
298 |
> |
sam_toks[i].find("XS:i:")==string::npos) |
299 |
> |
auxdata.push_back(sam_toks[i]); |
300 |
> |
|
301 |
> |
if (sam_toks[i].find("XS:A:")!=string::npos) |
302 |
> |
XS_found = true; |
303 |
> |
} |
304 |
> |
} |
305 |
> |
string aux("NH:i:"); |
306 |
> |
str_appendInt(aux, num_hits); |
307 |
> |
auxdata.push_back(aux); |
308 |
> |
if (next_hit) { |
309 |
> |
const char* nh_ref_name = "="; |
310 |
> |
nh_ref_name = rt.get_name(next_hit->ref_id()); |
311 |
> |
assert (nh_ref_name != NULL); |
312 |
> |
bool same_contig=(next_hit->ref_id()==bh.ref_id()); |
313 |
> |
aux="CC:Z:"; |
314 |
> |
aux+= (same_contig ? "=" : nh_ref_name); |
315 |
> |
auxdata.push_back(aux); |
316 |
> |
aux="CP:i:"; |
317 |
> |
int nh_gpos=next_hit->left() + 1; |
318 |
> |
str_appendInt(aux, nh_gpos); |
319 |
|
auxdata.push_back(aux); |
320 |
< |
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 |
320 |
> |
} //has next_hit |
321 |
|
// FIXME: this code is still a bit brittle, because it contains no |
322 |
|
// consistency check that the mates are on opposite strands - a current protocol |
323 |
|
// requirement, and that the strand indicated by the alignment is consistent |
324 |
|
// with the orientation of the splices (though that should be handled upstream). |
325 |
+ |
if (!XS_found) { |
326 |
|
const string xs_f("XS:A:+"); |
327 |
|
const string xs_r("XS:A:-"); |
328 |
< |
if (bh.contiguous()) { |
329 |
< |
if (library_type == FR_FIRSTSTRAND) { |
330 |
< |
if (insert_side == FRAG_LEFT || insert_side == FRAG_UNPAIRED) { |
331 |
< |
if (bh.antisense_align()) |
332 |
< |
auxdata.push_back(xs_f); |
333 |
< |
else |
334 |
< |
auxdata.push_back(xs_r); |
335 |
< |
} |
336 |
< |
else { |
337 |
< |
if (bh.antisense_align()) |
338 |
< |
auxdata.push_back(xs_r); |
339 |
< |
else |
340 |
< |
auxdata.push_back(xs_f); |
341 |
< |
} |
342 |
< |
} |
343 |
< |
else if (library_type == FR_SECONDSTRAND) { |
344 |
< |
if (insert_side == FRAG_LEFT || insert_side == FRAG_UNPAIRED){ |
345 |
< |
if (bh.antisense_align()) |
346 |
< |
auxdata.push_back(xs_r); |
347 |
< |
else |
348 |
< |
auxdata.push_back(xs_f); |
349 |
< |
} |
350 |
< |
else |
351 |
< |
{ |
352 |
< |
if (bh.antisense_align()) |
353 |
< |
auxdata.push_back(xs_f); |
354 |
< |
else |
355 |
< |
auxdata.push_back(xs_r); |
356 |
< |
} |
357 |
< |
} |
358 |
< |
} //bh.contiguous() |
243 |
< |
if (hitIndex >= 0) |
328 |
> |
if (library_type == FR_FIRSTSTRAND) { |
329 |
> |
if (insert_side == FRAG_LEFT || insert_side == FRAG_UNPAIRED) { |
330 |
> |
if (bh.antisense_align()) |
331 |
> |
auxdata.push_back(xs_f); |
332 |
> |
else |
333 |
> |
auxdata.push_back(xs_r); |
334 |
> |
} |
335 |
> |
else { |
336 |
> |
if (bh.antisense_align()) |
337 |
> |
auxdata.push_back(xs_r); |
338 |
> |
else |
339 |
> |
auxdata.push_back(xs_f); |
340 |
> |
} |
341 |
> |
} |
342 |
> |
else if (library_type == FR_SECONDSTRAND) { |
343 |
> |
if (insert_side == FRAG_LEFT || insert_side == FRAG_UNPAIRED){ |
344 |
> |
if (bh.antisense_align()) |
345 |
> |
auxdata.push_back(xs_r); |
346 |
> |
else |
347 |
> |
auxdata.push_back(xs_f); |
348 |
> |
} |
349 |
> |
else |
350 |
> |
{ |
351 |
> |
if (bh.antisense_align()) |
352 |
> |
auxdata.push_back(xs_f); |
353 |
> |
else |
354 |
> |
auxdata.push_back(xs_r); |
355 |
> |
} |
356 |
> |
} |
357 |
> |
} |
358 |
> |
if (hitIndex >= 0) |
359 |
|
{ |
360 |
< |
string aux("HI:i:"); |
361 |
< |
str_appendInt(aux, hitIndex); |
362 |
< |
auxdata.push_back(aux); |
360 |
> |
string aux("HI:i:"); |
361 |
> |
str_appendInt(aux, hitIndex); |
362 |
> |
auxdata.push_back(aux); |
363 |
|
} |
364 |
|
} |
365 |
|
|
368 |
|
const BowtieHit& bh, |
369 |
|
const char* bwt_buf, |
370 |
|
const char* read_alt_name, |
256 |
– |
const FragmentAlignmentGrade& grade, |
371 |
|
FragmentType insert_side, |
372 |
|
int num_hits, |
373 |
|
const BowtieHit* next_hit, |
380 |
|
tokenize(bwt_buf, "\t", sam_toks); |
381 |
|
|
382 |
|
string ref_name = sam_toks[2], ref_name2 = ""; |
383 |
< |
char rebuf2[2048] = {0}; |
270 |
< |
char cigar1[256] = {0}, cigar2[256] = {0}; |
383 |
> |
char cigar1[1024] = {0}, cigar2[1024] = {0}; |
384 |
|
string left_seq, right_seq, left_qual, right_qual; |
385 |
|
int left = -1, left1 = -1, left2 = -1; |
386 |
|
bool fusion_alignment = false; |
387 |
|
size_t XF_index = 0; |
388 |
|
for (size_t i = 11; i < sam_toks.size(); ++i) |
389 |
|
{ |
390 |
< |
const string& tok = sam_toks[i]; |
390 |
> |
string& tok = sam_toks[i]; |
391 |
|
if (strncmp(tok.c_str(), "XF", 2) == 0) |
392 |
|
{ |
393 |
|
fusion_alignment = true; |
409 |
|
|
410 |
|
break; |
411 |
|
} |
412 |
+ |
|
413 |
+ |
else if (strncmp(tok.c_str(), "AS", 2) == 0) |
414 |
+ |
{ |
415 |
+ |
char AS_score[128] = {0}; |
416 |
+ |
sprintf(AS_score, "AS:i:%d", bh.alignment_score()); |
417 |
+ |
tok = AS_score; |
418 |
+ |
} |
419 |
|
} |
420 |
|
|
421 |
|
string qname(read_alt_name); |
442 |
|
|
443 |
|
int gpos=isdigit(sam_toks[3][0]) ? atoi(sam_toks[3].c_str()) : 0; |
444 |
|
int mapQ=255; |
445 |
< |
if (grade.num_alignments > 1) { |
446 |
< |
double err_prob = 1 - (1.0 / grade.num_alignments); |
445 |
> |
if (num_hits > 1) { |
446 |
> |
double err_prob = 1 - (1.0 / num_hits); |
447 |
|
mapQ = (int)(-10.0 * log(err_prob) / log(10.0)); |
448 |
|
} |
449 |
|
int tlen =atoi(sam_toks[8].c_str()); //TLEN |
514 |
|
flag |= 0x100; |
515 |
|
|
516 |
|
string ref_name = sam_toks[2], ref_name2 = ""; |
517 |
< |
char rebuf2[2048] = {0}; |
398 |
< |
char cigar1[256] = {0}, cigar2[256] = {0}; |
517 |
> |
char cigar1[1024] = {0}, cigar2[1024] = {0}; |
518 |
|
string left_seq, right_seq, left_qual, right_qual; |
519 |
|
int left = -1, left1 = -1, left2 = -1; |
520 |
|
bool fusion_alignment = false; |
521 |
|
size_t XF_tok_idx = 11; |
522 |
|
for (; XF_tok_idx < sam_toks.size(); ++XF_tok_idx) |
523 |
|
{ |
524 |
< |
const string& tok = sam_toks[XF_tok_idx]; |
524 |
> |
string& tok = sam_toks[XF_tok_idx]; |
525 |
|
if (strncmp(tok.c_str(), "XF", 2) == 0) |
526 |
|
{ |
527 |
|
fusion_alignment = true; |
542 |
|
|
543 |
|
break; |
544 |
|
} |
545 |
+ |
|
546 |
+ |
else if (strncmp(tok.c_str(), "AS", 2) == 0) |
547 |
+ |
{ |
548 |
+ |
char AS_score[128] = {0}; |
549 |
+ |
sprintf(AS_score, "AS:i:%d", bh.alignment_score()); |
550 |
+ |
tok = AS_score; |
551 |
+ |
} |
552 |
|
} |
553 |
|
|
554 |
|
int gpos=isdigit(sam_toks[3][0]) ? atoi(sam_toks[3].c_str()) : 0; |
645 |
|
|
646 |
|
void print_sam_for_single(const RefSequenceTable& rt, |
647 |
|
const HitsForRead& hits, |
648 |
< |
const FragmentAlignmentGrade& grade, |
649 |
< |
FragmentType frag_type, |
650 |
< |
//FILE* reads_file, |
651 |
< |
Read& read, |
652 |
< |
GBamWriter& bam_writer, |
527 |
< |
FILE* um_out//write unmapped reads here |
528 |
< |
) |
648 |
> |
FragmentType frag_type, |
649 |
> |
Read& read, |
650 |
> |
GBamWriter& bam_writer, |
651 |
> |
GBamWriter* um_out //write unmapped reads here |
652 |
> |
) |
653 |
|
{ |
654 |
|
assert(!read.alt_name.empty()); |
655 |
+ |
if (hits.hits.empty()) |
656 |
+ |
return; |
657 |
+ |
|
658 |
|
lex_hit_sort s(rt, hits); |
659 |
|
vector<uint32_t> index_vector; |
660 |
|
|
664 |
|
|
665 |
|
sort(index_vector.begin(), index_vector.end(), s); |
666 |
|
|
667 |
< |
size_t primaryHit = (hits.hits.empty()? 0: random() % hits.hits.size()); |
667 |
> |
size_t primaryHit = 0; |
668 |
> |
if (!report_secondary_alignments) |
669 |
> |
primaryHit = random() % hits.hits.size(); |
670 |
> |
|
671 |
|
bool multipleHits = (hits.hits.size() > 1); |
672 |
|
for (i = 0; i < hits.hits.size(); ++i) |
673 |
|
{ |
544 |
– |
bool primary = (i == primaryHit); |
674 |
|
size_t index = index_vector[i]; |
675 |
+ |
bool primary = (index == primaryHit); |
676 |
|
const BowtieHit& bh = hits.hits[index]; |
677 |
|
rewrite_sam_record(bam_writer, rt, |
678 |
|
bh, |
679 |
|
bh.hitfile_rec().c_str(), |
680 |
|
read.alt_name.c_str(), |
551 |
– |
grade, |
681 |
|
frag_type, |
682 |
|
hits.hits.size(), |
683 |
|
(i < hits.hits.size()-1) ? &(hits.hits[index_vector[i+1]]) : NULL, |
687 |
|
} |
688 |
|
|
689 |
|
void print_sam_for_pair(const RefSequenceTable& rt, |
690 |
< |
const HitsForRead& left_hits, |
562 |
< |
const HitsForRead& right_hits, |
690 |
> |
const vector<pair<BowtieHit, BowtieHit> >& best_hits, |
691 |
|
const InsertAlignmentGrade& grade, |
692 |
|
ReadStream& left_reads_file, |
693 |
|
ReadStream& right_reads_file, |
694 |
|
GBamWriter& bam_writer, |
695 |
< |
FILE* left_um_out, |
696 |
< |
FILE* right_um_out, |
695 |
> |
GBamWriter* left_um_out, |
696 |
> |
GBamWriter* right_um_out, |
697 |
|
uint64_t begin_id = 0, |
698 |
|
uint64_t end_id = std::numeric_limits<uint64_t>::max()) |
699 |
|
{ |
572 |
– |
assert (left_hits.insert_id == right_hits.insert_id); |
573 |
– |
|
700 |
|
Read left_read; |
701 |
|
Read right_read; |
702 |
< |
assert (left_hits.hits.size() == right_hits.hits.size() || |
703 |
< |
(left_hits.hits.empty() || right_hits.hits.empty())); |
702 |
> |
if (best_hits.empty()) |
703 |
> |
return; |
704 |
> |
|
705 |
> |
size_t i; |
706 |
> |
HitsForRead right_hits; |
707 |
> |
for (i = 0; i < best_hits.size(); ++i) |
708 |
> |
right_hits.hits.push_back(best_hits[i].second); |
709 |
|
|
710 |
|
size_t primaryHit = 0; |
711 |
|
vector<uint32_t> index_vector; |
712 |
< |
if(right_hits.hits.size() > 0) |
713 |
< |
{ |
714 |
< |
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 |
< |
} |
712 |
> |
lex_hit_sort s(rt, right_hits); |
713 |
> |
for (i = 0; i < right_hits.hits.size(); ++i) |
714 |
> |
index_vector.push_back(i); |
715 |
|
|
716 |
< |
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 |
< |
} |
716 |
> |
sort(index_vector.begin(), index_vector.end(), s); |
717 |
|
|
718 |
< |
for (size_t i = 0; i < left_hits.hits.size(); ++i) |
719 |
< |
{ |
720 |
< |
bool primary = (i == primaryHit); |
721 |
< |
size_t index = index_vector[i]; |
722 |
< |
const BowtieHit& bh = left_hits.hits[index]; |
723 |
< |
rewrite_sam_record(bam_writer, rt, |
724 |
< |
bh, |
725 |
< |
bh.hitfile_rec().c_str(), |
726 |
< |
left_read.alt_name.c_str(), |
727 |
< |
grade, |
728 |
< |
FRAG_LEFT, |
729 |
< |
NULL, |
730 |
< |
left_hits.hits.size(), |
731 |
< |
(i < left_hits.hits.size() - 1) ? &(left_hits.hits[index_vector[i+1]]) : NULL, |
732 |
< |
primary, |
733 |
< |
-1); |
734 |
< |
|
735 |
< |
} |
736 |
< |
} |
737 |
< |
else |
738 |
< |
{ |
739 |
< |
assert (false); |
740 |
< |
} |
718 |
> |
if (!report_secondary_alignments) |
719 |
> |
primaryHit = random() % right_hits.hits.size(); |
720 |
> |
|
721 |
> |
bool got_left_read = left_reads_file.getRead(best_hits[0].first.insert_id(), left_read, |
722 |
> |
reads_format, false, begin_id, end_id, |
723 |
> |
left_um_out, false); |
724 |
> |
|
725 |
> |
bool got_right_read = right_reads_file.getRead(best_hits[0].first.insert_id(), right_read, |
726 |
> |
reads_format, false, begin_id, end_id, |
727 |
> |
right_um_out, false); |
728 |
> |
|
729 |
> |
assert (got_left_read && got_right_read); |
730 |
> |
bool multipleHits = (best_hits.size() > 1); |
731 |
> |
for (i = 0; i < best_hits.size(); ++i) |
732 |
> |
{ |
733 |
> |
size_t index = index_vector[i]; |
734 |
> |
bool primary = (index == primaryHit); |
735 |
> |
const BowtieHit& right_bh = best_hits[index].second; |
736 |
> |
const BowtieHit& left_bh = best_hits[index].first; |
737 |
> |
|
738 |
> |
rewrite_sam_record(bam_writer, rt, |
739 |
> |
right_bh, |
740 |
> |
right_bh.hitfile_rec().c_str(), |
741 |
> |
right_read.alt_name.c_str(), |
742 |
> |
grade, |
743 |
> |
FRAG_RIGHT, |
744 |
> |
&left_bh, |
745 |
> |
best_hits.size(), |
746 |
> |
(i < best_hits.size() - 1) ? &(best_hits[index_vector[i+1]].second) : NULL, |
747 |
> |
primary, |
748 |
> |
(multipleHits? i: -1)); |
749 |
> |
rewrite_sam_record(bam_writer, rt, |
750 |
> |
left_bh, |
751 |
> |
left_bh.hitfile_rec().c_str(), |
752 |
> |
left_read.alt_name.c_str(), |
753 |
> |
grade, |
754 |
> |
FRAG_LEFT, |
755 |
> |
&right_bh, |
756 |
> |
best_hits.size(), |
757 |
> |
(i < best_hits.size() - 1) ? &(best_hits[index_vector[i+1]].first) : NULL, |
758 |
> |
primary, |
759 |
> |
(multipleHits? i: -1)); |
760 |
> |
} |
761 |
|
} |
762 |
|
|
763 |
|
/** |
770 |
|
InsertionSet& insertions, |
771 |
|
DeletionSet& deletions) |
772 |
|
{ |
773 |
< |
for (size_t i = 0; i < hits.hits.size(); ++i) |
773 |
> |
for (size_t i = 0; i < hits.hits.size(); ++i) |
774 |
> |
{ |
775 |
> |
const BowtieHit& bh = hits.hits[i]; |
776 |
> |
insertions_from_alignment(bh, insertions); |
777 |
> |
deletions_from_alignment(bh, deletions); |
778 |
> |
} |
779 |
> |
} |
780 |
> |
|
781 |
> |
void update_coverage(const HitsForRead& hits, |
782 |
> |
Coverage& coverage) |
783 |
> |
{ |
784 |
> |
for (size_t i = 0; i < hits.hits.size(); ++i) |
785 |
> |
{ |
786 |
> |
const BowtieHit& hit = hits.hits[i]; |
787 |
> |
const vector<CigarOp>& cigar = hit.cigar(); |
788 |
> |
unsigned int positionInGenome = hit.left(); |
789 |
> |
RefID ref_id = hit.ref_id(); |
790 |
> |
|
791 |
> |
for(size_t c = 0; c < cigar.size(); ++c) |
792 |
|
{ |
793 |
< |
const BowtieHit& bh = hits.hits[i]; |
794 |
< |
insertions_from_alignment(bh, insertions); |
795 |
< |
deletions_from_alignment(bh, deletions); |
793 |
> |
int opcode = cigar[c].opcode; |
794 |
> |
int length = cigar[c].length; |
795 |
> |
switch(opcode) |
796 |
> |
{ |
797 |
> |
case REF_SKIP: |
798 |
> |
case MATCH: |
799 |
> |
case DEL: |
800 |
> |
if (opcode == MATCH) |
801 |
> |
coverage.add_coverage(ref_id, positionInGenome, length); |
802 |
> |
|
803 |
> |
positionInGenome += length; |
804 |
> |
break; |
805 |
> |
case rEF_SKIP: |
806 |
> |
case mATCH: |
807 |
> |
case dEL: |
808 |
> |
positionInGenome -= length; |
809 |
> |
|
810 |
> |
if (opcode == mATCH) |
811 |
> |
coverage.add_coverage(ref_id, positionInGenome + 1, length); |
812 |
> |
break; |
813 |
> |
case FUSION_FF: |
814 |
> |
case FUSION_FR: |
815 |
> |
case FUSION_RF: |
816 |
> |
case FUSION_RR: |
817 |
> |
positionInGenome = length; |
818 |
> |
ref_id = hit.ref_id2(); |
819 |
> |
break; |
820 |
> |
default: |
821 |
> |
break; |
822 |
> |
} |
823 |
|
} |
824 |
+ |
} |
825 |
|
} |
826 |
|
|
827 |
|
|
722 |
– |
FusionSet empty_fusions; |
828 |
|
void update_fusions(const HitsForRead& hits, |
829 |
|
RefSequenceTable& rt, |
830 |
|
FusionSet& fusions, |
831 |
< |
FusionSet& fusions_ref = empty_fusions) |
831 |
> |
const FusionSet& fusions_ref = empty_fusions) |
832 |
|
{ |
833 |
+ |
if (!fusion_search) |
834 |
+ |
return; |
835 |
+ |
|
836 |
|
if (hits.hits.size() > fusion_multireads) |
837 |
|
return; |
838 |
|
|
844 |
|
if (bh.edit_dist() > fusion_read_mismatches) |
845 |
|
continue; |
846 |
|
|
739 |
– |
// daehwan |
740 |
– |
#if 0 |
741 |
– |
vector<Fusion> new_fusions; |
742 |
– |
fusions_from_spliced_hit(bh, new_fusions); |
743 |
– |
|
744 |
– |
for(size_t i = 0; i < new_fusions.size(); ++i) |
745 |
– |
{ |
746 |
– |
Fusion fusion = new_fusions[i]; |
747 |
– |
if (fusion.left == 110343414 && fusion.right == 80211173 && hits.hits.size() > 1) |
748 |
– |
{ |
749 |
– |
for (size_t t = 0; t < hits.hits.size(); ++t) |
750 |
– |
{ |
751 |
– |
const BowtieHit& lh = hits.hits[t]; |
752 |
– |
cout << "insert id: " << lh.insert_id() << endl; |
753 |
– |
cout << (lh.antisense_align() ? "-" : "+") << endl; |
754 |
– |
cout << lh.ref_id() << ": " << lh.left() << "-" << lh.right() << endl; |
755 |
– |
cout << lh.ref_id2() << ": " << print_cigar(lh.cigar()) << endl; |
756 |
– |
} |
757 |
– |
} |
758 |
– |
} |
759 |
– |
#endif |
760 |
– |
|
847 |
|
fusions_from_alignment(bh, fusions, rt, update_stat); |
848 |
|
|
849 |
|
if (update_stat) |
854 |
|
void update_junctions(const HitsForRead& hits, |
855 |
|
JunctionSet& junctions) |
856 |
|
{ |
857 |
< |
for (size_t i = 0; i < hits.hits.size(); ++i) |
858 |
< |
{ |
859 |
< |
const BowtieHit& bh = hits.hits[i]; |
860 |
< |
junctions_from_alignment(bh, junctions); |
861 |
< |
} |
857 |
> |
for (size_t i = 0; i < hits.hits.size(); ++i) |
858 |
> |
{ |
859 |
> |
const BowtieHit& bh = hits.hits[i]; |
860 |
> |
junctions_from_alignment(bh, junctions); |
861 |
> |
} |
862 |
|
} |
863 |
|
|
864 |
|
// Extracts junctions from all the SAM hits (based on REF_SKIPs) in the hit file |
871 |
|
for (size_t i = 0; i < hits.hits.size(); ++i) |
872 |
|
{ |
873 |
|
BowtieHit& bh = hits.hits[i]; |
874 |
+ |
if (bh.edit_dist() > max_read_mismatches) |
875 |
+ |
continue; |
876 |
+ |
|
877 |
|
bool filter_hit = false; |
878 |
|
if (!bh.contiguous()) |
879 |
|
{ |
935 |
|
{ |
936 |
|
HitsForRead best_hits; |
937 |
|
best_hits.insert_id = curr_left_obs_order; |
849 |
– |
FragmentAlignmentGrade grade; |
938 |
|
|
939 |
|
// Process hits for left singleton, select best alignments |
940 |
< |
read_best_alignments(curr_left_hit_group, grade, best_hits, *gtf_junctions); |
940 |
> |
read_best_alignments(curr_left_hit_group, best_hits, *gtf_junctions); |
941 |
> |
update_coverage(best_hits, *coverage); |
942 |
|
update_junctions(best_hits, *junctions); |
943 |
+ |
update_insertions_and_deletions(best_hits, *insertions, *deletions); |
944 |
|
update_fusions(best_hits, *rt, *fusions); |
945 |
< |
|
945 |
> |
|
946 |
|
// Get next hit group |
947 |
|
l_hs.next_read_hits(curr_left_hit_group); |
948 |
|
curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id); |
957 |
|
best_hits.insert_id = curr_right_obs_order; |
958 |
|
if (curr_right_obs_order >= begin_id) |
959 |
|
{ |
870 |
– |
FragmentAlignmentGrade grade; |
960 |
|
// Process hit for right singleton, select best alignments |
961 |
< |
read_best_alignments(curr_right_hit_group,grade, best_hits, *gtf_junctions); |
961 |
> |
read_best_alignments(curr_right_hit_group, best_hits, *gtf_junctions); |
962 |
> |
update_coverage(best_hits, *coverage); |
963 |
|
update_junctions(best_hits, *junctions); |
964 |
+ |
update_insertions_and_deletions(best_hits, *insertions, *deletions); |
965 |
|
update_fusions(best_hits, *rt, *fusions); |
966 |
|
} |
967 |
|
|
976 |
|
curr_left_obs_order < end_id && |
977 |
|
curr_left_obs_order != VMAXINT32) |
978 |
|
{ |
979 |
< |
if (curr_left_hit_group.hits.empty()) |
980 |
< |
{ |
981 |
< |
HitsForRead right_best_hits; |
982 |
< |
right_best_hits.insert_id = curr_right_obs_order; |
983 |
< |
|
984 |
< |
FragmentAlignmentGrade grade; |
985 |
< |
read_best_alignments(curr_right_hit_group, grade, right_best_hits, *gtf_junctions); |
986 |
< |
update_junctions(right_best_hits, *junctions); |
987 |
< |
update_fusions(right_best_hits, *rt, *fusions); |
988 |
< |
} |
989 |
< |
else if (curr_right_hit_group.hits.empty()) |
979 |
> |
vector<pair<BowtieHit, BowtieHit> > best_hits; |
980 |
> |
|
981 |
> |
InsertAlignmentGrade grade; |
982 |
> |
pair_best_alignments(curr_left_hit_group, |
983 |
> |
curr_right_hit_group, |
984 |
> |
grade, |
985 |
> |
best_hits, |
986 |
> |
*gtf_junctions); |
987 |
> |
|
988 |
> |
HitsForRead best_left_hit_group; best_left_hit_group.insert_id = curr_left_obs_order; |
989 |
> |
HitsForRead best_right_hit_group; best_right_hit_group.insert_id = curr_left_obs_order; |
990 |
> |
|
991 |
> |
for (size_t i = 0; i < best_hits.size(); ++i) |
992 |
|
{ |
993 |
< |
HitsForRead left_best_hits; |
994 |
< |
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); |
993 |
> |
best_left_hit_group.hits.push_back(best_hits[i].first); |
994 |
> |
best_right_hit_group.hits.push_back(best_hits[i].second); |
995 |
|
} |
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! |
996 |
|
|
997 |
< |
InsertAlignmentGrade grade; |
998 |
< |
pair_best_alignments(curr_left_hit_group, |
999 |
< |
curr_right_hit_group, |
1000 |
< |
grade, |
1001 |
< |
left_best_hits, |
1002 |
< |
right_best_hits); |
1003 |
< |
update_junctions(left_best_hits, *junctions); |
1004 |
< |
update_junctions(right_best_hits, *junctions); |
1005 |
< |
update_fusions(left_best_hits, *rt, *fusions); |
1006 |
< |
update_fusions(right_best_hits, *rt, *fusions); |
928 |
< |
} |
929 |
< |
|
997 |
> |
update_coverage(best_left_hit_group, *coverage); |
998 |
> |
update_junctions(best_left_hit_group, *junctions); |
999 |
> |
update_insertions_and_deletions(best_left_hit_group, *insertions, *deletions); |
1000 |
> |
update_fusions(best_left_hit_group, *rt, *fusions); |
1001 |
> |
|
1002 |
> |
update_coverage(best_right_hit_group, *coverage); |
1003 |
> |
update_junctions(best_right_hit_group, *junctions); |
1004 |
> |
update_insertions_and_deletions(best_right_hit_group, *insertions, *deletions); |
1005 |
> |
update_fusions(best_right_hit_group, *rt, *fusions); |
1006 |
> |
|
1007 |
|
l_hs.next_read_hits(curr_left_hit_group); |
1008 |
|
curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id); |
1009 |
|
|
1021 |
|
string left_map_fname; |
1022 |
|
string right_map_fname; |
1023 |
|
|
947 |
– |
JunctionSet* junctions; |
948 |
– |
JunctionSet* gtf_junctions; |
949 |
– |
FusionSet* fusions; |
1024 |
|
RefSequenceTable* rt; |
1025 |
|
|
1026 |
+ |
JunctionSet* gtf_junctions; |
1027 |
+ |
|
1028 |
|
uint64_t begin_id; |
1029 |
|
uint64_t end_id; |
1030 |
|
int64_t left_map_offset; |
1031 |
|
int64_t right_map_offset; |
1032 |
+ |
|
1033 |
+ |
JunctionSet* junctions; |
1034 |
+ |
InsertionSet* insertions; |
1035 |
+ |
DeletionSet* deletions; |
1036 |
+ |
FusionSet* fusions; |
1037 |
+ |
|
1038 |
+ |
Coverage* coverage; |
1039 |
|
}; |
1040 |
|
|
1041 |
|
struct ReportWorker |
1048 |
|
ReadStream left_reads_file(left_reads_fname); |
1049 |
|
if (left_reads_file.file() == NULL) |
1050 |
|
err_die("Error: cannot open %s for reading\n", left_reads_fname.c_str()); |
1051 |
< |
|
1051 |
> |
if (left_reads_file.isBam()) { |
1052 |
> |
left_reads_file.use_alt_name(); |
1053 |
> |
left_reads_file.ignoreQC(); |
1054 |
> |
} |
1055 |
|
if (left_reads_offset > 0) |
1056 |
|
left_reads_file.seek(left_reads_offset); |
1057 |
|
|
1058 |
< |
if (!zpacker.empty()) left_um_fname+=".z"; |
1059 |
< |
FZPipe left_um_out; |
1060 |
< |
if (left_um_out.openWrite(left_um_fname.c_str(), zpacker)==NULL) |
975 |
< |
err_die("Error: cannot open file %s for writing!\n",left_um_fname.c_str()); |
1058 |
> |
//if (!zpacker.empty()) left_um_fname+=".z"; |
1059 |
> |
GBamWriter* left_um_out=new GBamWriter(left_um_fname.c_str(), sam_header.c_str()); |
1060 |
> |
GBamWriter* right_um_out=NULL; |
1061 |
|
|
1062 |
+ |
//if (left_um_out.openWrite(left_um_fname.c_str(), zpacker)==NULL) |
1063 |
+ |
// err_die("Error: cannot open file %s for writing!\n",left_um_fname.c_str()); |
1064 |
|
ReadStream right_reads_file(right_reads_fname); |
1065 |
|
if (right_reads_offset > 0) |
1066 |
|
right_reads_file.seek(right_reads_offset); |
1067 |
|
|
1068 |
< |
FZPipe right_um_out; |
1069 |
< |
if (right_reads_fname != "") |
1068 |
> |
//FZPipe right_um_out; |
1069 |
> |
if (!right_reads_fname.empty()) |
1070 |
|
{ |
1071 |
< |
if (!zpacker.empty()) right_um_fname+=".z"; |
1072 |
< |
if (right_um_out.openWrite(right_um_fname.c_str(), zpacker)==NULL) |
1073 |
< |
err_die("Error: cannot open file %s for writing!\n",right_um_fname.c_str()); |
1071 |
> |
if (right_reads_file.isBam()) { |
1072 |
> |
right_reads_file.use_alt_name(); |
1073 |
> |
right_reads_file.ignoreQC(); |
1074 |
> |
right_um_out=new GBamWriter(right_um_fname.c_str(), sam_header.c_str()); |
1075 |
> |
} |
1076 |
> |
//if (!zpacker.empty()) right_um_fname+=".z"; |
1077 |
> |
//if (right_um_out.openWrite(right_um_fname.c_str(), zpacker)==NULL) |
1078 |
> |
// err_die("Error: cannot open file %s for writing!\n",right_um_fname.c_str()); |
1079 |
|
} |
1080 |
|
|
1081 |
|
vector<BAMHitFactory*> hit_factories; |
1098 |
|
uint64_t curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id); |
1099 |
|
uint64_t curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id); |
1100 |
|
|
1101 |
+ |
const bool final_report = true; |
1102 |
+ |
|
1103 |
|
// While we still have unreported hits... |
1104 |
|
Read l_read; |
1105 |
|
Read r_read; |
1113 |
|
{ |
1114 |
|
HitsForRead best_hits; |
1115 |
|
best_hits.insert_id = curr_left_obs_order; |
1022 |
– |
FragmentAlignmentGrade grade; |
1116 |
|
bool got_read=left_reads_file.getRead(curr_left_obs_order, l_read, |
1117 |
< |
reads_format, false, left_um_out.file, |
1118 |
< |
false, begin_id, end_id); |
1119 |
< |
assert(got_read); |
1117 |
> |
reads_format, false, begin_id, end_id, |
1118 |
> |
left_um_out, false); |
1119 |
> |
//assert(got_read); |
1120 |
|
|
1121 |
|
if (right_reads_file.file()) { |
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()); |
1122 |
|
got_read=right_reads_file.getRead(curr_left_obs_order, r_read, |
1123 |
< |
reads_format, false, right_um_out.file, |
1124 |
< |
true, begin_id, end_id); |
1125 |
< |
assert(got_read); |
1123 |
> |
reads_format, false, begin_id, end_id, |
1124 |
> |
right_um_out, true); |
1125 |
> |
//assert(got_read); |
1126 |
|
} |
1127 |
< |
|
1127 |
> |
|
1128 |
|
exclude_hits_on_filtered_junctions(*junctions, curr_left_hit_group); |
1129 |
|
|
1130 |
|
// Process hits for left singleton, select best alignments |
1131 |
< |
read_best_alignments(curr_left_hit_group, grade, best_hits, *gtf_junctions); |
1131 |
> |
read_best_alignments(curr_left_hit_group, best_hits, *gtf_junctions, |
1132 |
> |
*junctions, *insertions, *deletions, *fusions, *coverage, |
1133 |
> |
final_report); |
1134 |
|
|
1135 |
< |
if (best_hits.hits.size()>0 && best_hits.hits.size() <= max_multihits) |
1135 |
> |
if (best_hits.hits.size() > 0) |
1136 |
|
{ |
1137 |
|
update_junctions(best_hits, *final_junctions); |
1138 |
|
update_insertions_and_deletions(best_hits, *final_insertions, *final_deletions); |
1140 |
|
|
1141 |
|
print_sam_for_single(*rt, |
1142 |
|
best_hits, |
1050 |
– |
grade, |
1143 |
|
(right_map_fname.empty() ? FRAG_UNPAIRED : FRAG_LEFT), |
1144 |
|
l_read, |
1145 |
|
bam_writer, |
1146 |
< |
left_um_out.file); |
1146 |
> |
left_um_out); |
1147 |
|
} |
1148 |
|
// Get next hit group |
1149 |
|
left_hs.next_read_hits(curr_left_hit_group); |
1157 |
|
{ |
1158 |
|
HitsForRead best_hits; |
1159 |
|
best_hits.insert_id = curr_right_obs_order; |
1068 |
– |
FragmentAlignmentGrade grade; |
1160 |
|
|
1161 |
|
if (curr_right_obs_order >= begin_id) |
1162 |
|
{ |
1163 |
|
bool got_read=right_reads_file.getRead(curr_right_obs_order, r_read, |
1164 |
< |
reads_format, false, right_um_out.file, |
1165 |
< |
false, begin_id, end_id); |
1075 |
< |
|
1076 |
< |
assert(got_read); |
1164 |
> |
reads_format, false, begin_id, end_id, |
1165 |
> |
right_um_out, false); |
1166 |
|
|
1167 |
+ |
//assert(got_read); |
1168 |
+ |
/* |
1169 |
|
fprintf(right_um_out.file, "@%s #MAPPED#\n%s\n+\n%s\n", r_read.alt_name.c_str(), |
1170 |
|
r_read.seq.c_str(), r_read.qual.c_str()); |
1171 |
+ |
*/ |
1172 |
|
got_read=left_reads_file.getRead(curr_right_obs_order, l_read, |
1173 |
< |
reads_format, false, left_um_out.file, |
1174 |
< |
true, begin_id, end_id); |
1175 |
< |
assert(got_read); |
1173 |
> |
reads_format, false, begin_id, end_id, |
1174 |
> |
left_um_out, true); |
1175 |
> |
//assert(got_read); |
1176 |
|
|
1177 |
|
exclude_hits_on_filtered_junctions(*junctions, curr_right_hit_group); |
1178 |
|
|
1179 |
|
// Process hit for right singleton, select best alignments |
1180 |
< |
read_best_alignments(curr_right_hit_group, grade, best_hits, *gtf_junctions); |
1181 |
< |
if (best_hits.hits.size()>0 && best_hits.hits.size() <= max_multihits) |
1180 |
> |
read_best_alignments(curr_right_hit_group, best_hits, *gtf_junctions, |
1181 |
> |
*junctions, *insertions, *deletions, *fusions, *coverage, |
1182 |
> |
final_report); |
1183 |
> |
|
1184 |
> |
if (best_hits.hits.size() > 0) |
1185 |
|
{ |
1186 |
|
update_junctions(best_hits, *final_junctions); |
1187 |
|
update_insertions_and_deletions(best_hits, *final_insertions, *final_deletions); |
1189 |
|
|
1190 |
|
print_sam_for_single(*rt, |
1191 |
|
best_hits, |
1097 |
– |
grade, |
1192 |
|
FRAG_RIGHT, |
1193 |
|
r_read, |
1194 |
|
bam_writer, |
1195 |
< |
right_um_out.file); |
1195 |
> |
right_um_out); |
1196 |
|
} |
1197 |
|
} |
1198 |
|
|
1214 |
|
//write it in the mapped file with the #MAPPED# flag |
1215 |
|
|
1216 |
|
bool got_read=right_reads_file.getRead(curr_left_obs_order, r_read, |
1217 |
< |
reads_format, false, right_um_out.file, |
1218 |
< |
false, begin_id, end_id); |
1219 |
< |
assert(got_read); |
1217 |
> |
reads_format, false, begin_id, end_id, |
1218 |
> |
right_um_out, false); |
1219 |
> |
//assert(got_read); |
1220 |
> |
/* |
1221 |
|
fprintf(right_um_out.file, "@%s #MAPPED#\n%s\n+\n%s\n", r_read.alt_name.c_str(), |
1222 |
|
r_read.seq.c_str(), r_read.qual.c_str()); |
1223 |
+ |
*/ |
1224 |
|
HitsForRead right_best_hits; |
1225 |
|
right_best_hits.insert_id = curr_right_obs_order; |
1226 |
|
|
1227 |
< |
FragmentAlignmentGrade grade; |
1228 |
< |
read_best_alignments(curr_right_hit_group, grade, right_best_hits, *gtf_junctions); |
1229 |
< |
|
1230 |
< |
if (right_best_hits.hits.size()>0 && right_best_hits.hits.size() <= max_multihits) |
1227 |
> |
read_best_alignments(curr_right_hit_group, right_best_hits, *gtf_junctions, |
1228 |
> |
*junctions, *insertions, *deletions, *fusions, *coverage, |
1229 |
> |
final_report); |
1230 |
> |
|
1231 |
> |
if (right_best_hits.hits.size() > 0) |
1232 |
|
{ |
1233 |
|
update_junctions(right_best_hits, *final_junctions); |
1234 |
|
update_insertions_and_deletions(right_best_hits, *final_insertions, *final_deletions); |
1236 |
|
|
1237 |
|
print_sam_for_single(*rt, |
1238 |
|
right_best_hits, |
1142 |
– |
grade, |
1239 |
|
FRAG_RIGHT, |
1240 |
|
r_read, |
1241 |
|
bam_writer, |
1242 |
< |
right_um_out.file); |
1242 |
> |
right_um_out); |
1243 |
|
} |
1244 |
|
} |
1245 |
|
else if (curr_right_hit_group.hits.empty()) |
1248 |
|
left_best_hits.insert_id = curr_left_obs_order; |
1249 |
|
//only left read mapped |
1250 |
|
bool got_read=left_reads_file.getRead(curr_left_obs_order, l_read, |
1251 |
< |
reads_format, false, left_um_out.file, |
1252 |
< |
false, begin_id, end_id); |
1253 |
< |
assert(got_read); |
1158 |
< |
fprintf(left_um_out.file, "@%s #MAPPED#\n%s\n+\n%s\n", l_read.alt_name.c_str(), |
1159 |
< |
l_read.seq.c_str(), l_read.qual.c_str()); |
1160 |
< |
FragmentAlignmentGrade grade; |
1251 |
> |
reads_format, false, begin_id, end_id, |
1252 |
> |
left_um_out, false); |
1253 |
> |
|
1254 |
|
// Process hits for left singleton, select best alignments |
1255 |
< |
read_best_alignments(curr_left_hit_group, grade, left_best_hits, *gtf_junctions); |
1256 |
< |
if (left_best_hits.hits.size()>0 && left_best_hits.hits.size() <= max_multihits) |
1255 |
> |
read_best_alignments(curr_left_hit_group, left_best_hits, *gtf_junctions, |
1256 |
> |
*junctions, *insertions, *deletions, *fusions, *coverage, |
1257 |
> |
final_report); |
1258 |
> |
|
1259 |
> |
if (left_best_hits.hits.size() > 0) |
1260 |
|
{ |
1261 |
|
update_junctions(left_best_hits, *final_junctions); |
1262 |
|
update_insertions_and_deletions(left_best_hits, *final_insertions, *final_deletions); |
1264 |
|
|
1265 |
|
print_sam_for_single(*rt, |
1266 |
|
left_best_hits, |
1171 |
– |
grade, |
1267 |
|
FRAG_LEFT, |
1268 |
|
l_read, |
1269 |
|
bam_writer, |
1270 |
< |
left_um_out.file); |
1270 |
> |
left_um_out); |
1271 |
|
} |
1272 |
|
} |
1273 |
|
else |
1274 |
< |
{ //hits for both left and right reads |
1275 |
< |
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; |
1274 |
> |
{ |
1275 |
> |
vector<pair<BowtieHit, BowtieHit> > best_hits; |
1276 |
|
|
1277 |
|
InsertAlignmentGrade grade; |
1278 |
|
pair_best_alignments(curr_left_hit_group, |
1279 |
|
curr_right_hit_group, |
1280 |
< |
grade, |
1281 |
< |
left_best_hits, |
1282 |
< |
right_best_hits); |
1280 |
> |
grade, best_hits, *gtf_junctions, |
1281 |
> |
*junctions, *insertions, *deletions, *fusions, *coverage, |
1282 |
> |
final_report); |
1283 |
> |
|
1284 |
> |
HitsForRead best_left_hit_group; best_left_hit_group.insert_id = curr_left_obs_order; |
1285 |
> |
HitsForRead best_right_hit_group; best_right_hit_group.insert_id = curr_left_obs_order; |
1286 |
|
|
1287 |
< |
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) |
1287 |
> |
for (size_t i = 0; i < best_hits.size(); ++i) |
1288 |
|
{ |
1289 |
< |
update_junctions(left_best_hits, *final_junctions); |
1290 |
< |
update_junctions(right_best_hits, *final_junctions); |
1291 |
< |
update_insertions_and_deletions(left_best_hits, *final_insertions, *final_deletions); |
1198 |
< |
update_insertions_and_deletions(right_best_hits, *final_insertions, *final_deletions); |
1289 |
> |
best_left_hit_group.hits.push_back(best_hits[i].first); |
1290 |
> |
best_right_hit_group.hits.push_back(best_hits[i].second); |
1291 |
> |
} |
1292 |
|
|
1293 |
< |
pair_support(left_best_hits, right_best_hits, *final_fusions, *fusions); |
1294 |
< |
update_fusions(left_best_hits, *rt, *final_fusions, *fusions); |
1295 |
< |
update_fusions(right_best_hits, *rt, *final_fusions, *fusions); |
1293 |
> |
if (best_hits.size() > 0) |
1294 |
> |
{ |
1295 |
> |
update_junctions(best_left_hit_group, *final_junctions); |
1296 |
> |
update_insertions_and_deletions(best_left_hit_group, *final_insertions, *final_deletions); |
1297 |
> |
update_fusions(best_left_hit_group, *rt, *final_fusions, *fusions); |
1298 |
> |
|
1299 |
> |
update_junctions(best_right_hit_group, *final_junctions); |
1300 |
> |
update_insertions_and_deletions(best_right_hit_group, *final_insertions, *final_deletions); |
1301 |
> |
update_fusions(best_right_hit_group, *rt, *final_fusions, *fusions); |
1302 |
> |
|
1303 |
> |
pair_support(best_hits, *final_fusions, *fusions); |
1304 |
|
|
1305 |
|
print_sam_for_pair(*rt, |
1306 |
< |
left_best_hits, |
1206 |
< |
right_best_hits, |
1306 |
> |
best_hits, |
1307 |
|
grade, |
1308 |
|
left_reads_file, |
1309 |
|
right_reads_file, |
1310 |
|
bam_writer, |
1311 |
< |
left_um_out.file, |
1312 |
< |
right_um_out.file, |
1311 |
> |
left_um_out, |
1312 |
> |
right_um_out, |
1313 |
|
begin_id, |
1314 |
|
end_id); |
1315 |
|
} |
1323 |
|
} |
1324 |
|
} //while we still have unreported hits.. |
1325 |
|
//print the remaining unmapped reads at the end of each reads' stream |
1326 |
+ |
|
1327 |
|
left_reads_file.getRead(VMAXINT32, l_read, |
1328 |
< |
reads_format, false, left_um_out.file, |
1329 |
< |
false, begin_id, end_id); |
1328 |
> |
reads_format, false, begin_id, end_id, |
1329 |
> |
left_um_out, false); |
1330 |
|
if (right_reads_file.file()) |
1331 |
|
right_reads_file.getRead(VMAXINT32, r_read, |
1332 |
< |
reads_format, false, right_um_out.file, |
1333 |
< |
false, begin_id, end_id); |
1332 |
> |
reads_format, false, begin_id, end_id, |
1333 |
> |
right_um_out, false); |
1334 |
> |
|
1335 |
> |
|
1336 |
> |
// pclose (pipe close), which waits for a process to end, seems to conflict with boost::thread::join somehow, |
1337 |
> |
// resulting in deadlock like behavior. |
1338 |
> |
delete left_um_out; |
1339 |
> |
delete right_um_out; |
1340 |
|
|
1234 |
– |
left_um_out.close(); |
1235 |
– |
right_um_out.close(); |
1341 |
|
|
1342 |
|
for (size_t i = 0; i < hit_factories.size(); ++i) |
1343 |
|
delete hit_factories[i]; |
1357 |
|
string right_um_fname; |
1358 |
|
|
1359 |
|
JunctionSet* gtf_junctions; |
1360 |
+ |
|
1361 |
|
JunctionSet* junctions; |
1362 |
+ |
InsertionSet* insertions; |
1363 |
+ |
DeletionSet* deletions; |
1364 |
|
FusionSet* fusions; |
1365 |
+ |
Coverage* coverage; |
1366 |
+ |
|
1367 |
|
JunctionSet* final_junctions; |
1368 |
|
InsertionSet* final_insertions; |
1369 |
|
DeletionSet* final_deletions; |
1403 |
|
JunctionSet gtf_junctions; |
1404 |
|
if (!gtf_juncs.empty()) |
1405 |
|
{ |
1406 |
< |
char splice_buf[2048]; |
1406 |
> |
char splice_buf[4096]; |
1407 |
|
FILE* splice_coords = fopen(gtf_juncs.c_str(), "r"); |
1408 |
|
if (splice_coords) |
1409 |
|
{ |
1430 |
|
uint32_t right_coord = atoi(scan_right_coord); |
1431 |
|
bool antisense = *orientation == '-'; |
1432 |
|
|
1433 |
< |
// add 1 to left_coord to meet BED format |
1324 |
< |
gtf_junctions.insert(make_pair<Junction, JunctionStats>(Junction(ref_id, left_coord + 1, right_coord, antisense), JunctionStats())); |
1433 |
> |
gtf_junctions.insert(make_pair<Junction, JunctionStats>(Junction(ref_id, left_coord, right_coord, antisense), JunctionStats())); |
1434 |
|
} |
1435 |
|
} |
1436 |
|
fprintf(stderr, "Loaded %d GFF junctions from %s.\n", (int)(gtf_junctions.size()), gtf_juncs.c_str()); |
1454 |
|
} |
1455 |
|
|
1456 |
|
JunctionSet vjunctions[num_threads]; |
1457 |
+ |
InsertionSet vinsertions[num_threads]; |
1458 |
+ |
DeletionSet vdeletions[num_threads]; |
1459 |
|
FusionSet vfusions[num_threads]; |
1460 |
+ |
Coverage vcoverages[num_threads]; |
1461 |
+ |
|
1462 |
|
vector<boost::thread*> threads; |
1463 |
|
for (int i = 0; i < num_threads; ++i) |
1464 |
|
{ |
1466 |
|
|
1467 |
|
worker.left_map_fname = left_map_fname; |
1468 |
|
worker.right_map_fname = right_map_fname; |
1469 |
+ |
worker.rt = &rt; |
1470 |
+ |
worker.gtf_junctions = >f_junctions; |
1471 |
+ |
|
1472 |
|
worker.junctions = &vjunctions[i]; |
1473 |
+ |
worker.insertions = &vinsertions[i]; |
1474 |
+ |
worker.deletions = &vdeletions[i]; |
1475 |
|
worker.fusions = &vfusions[i]; |
1476 |
< |
worker.gtf_junctions = >f_junctions; |
1359 |
< |
worker.rt = &rt; |
1476 |
> |
worker.coverage = &vcoverages[i]; |
1477 |
|
|
1478 |
|
worker.right_map_offset = 0; |
1479 |
|
|
1508 |
|
} |
1509 |
|
threads.clear(); |
1510 |
|
|
1511 |
< |
JunctionSet junctions = vjunctions[0]; |
1512 |
< |
FusionSet fusions = vfusions[0]; |
1511 |
> |
JunctionSet& junctions = vjunctions[0]; |
1512 |
> |
InsertionSet& insertions = vinsertions[0]; |
1513 |
> |
DeletionSet& deletions = vdeletions[0]; |
1514 |
> |
FusionSet& fusions = vfusions[0]; |
1515 |
> |
Coverage& coverage = vcoverages[0]; |
1516 |
|
for (size_t i = 1; i < num_threads; ++i) |
1517 |
|
{ |
1518 |
|
merge_with(junctions, vjunctions[i]); |
1519 |
|
vjunctions[i].clear(); |
1520 |
|
|
1521 |
+ |
merge_with(insertions, vinsertions[i]); |
1522 |
+ |
vinsertions[i].clear(); |
1523 |
+ |
|
1524 |
+ |
merge_with(deletions, vdeletions[i]); |
1525 |
+ |
vdeletions[i].clear(); |
1526 |
+ |
|
1527 |
|
merge_with(fusions, vfusions[i]); |
1528 |
|
vfusions[i].clear(); |
1529 |
+ |
|
1530 |
+ |
coverage.merge_with(vcoverages[i]); |
1531 |
+ |
vcoverages[i].clear(); |
1532 |
|
} |
1533 |
|
|
1534 |
+ |
coverage.calculate_coverage(); |
1535 |
+ |
|
1536 |
|
size_t num_unfiltered_juncs = junctions.size(); |
1537 |
|
fprintf(stderr, "Loaded %lu junctions\n", (long unsigned int) num_unfiltered_juncs); |
1538 |
|
|
1539 |
|
// Read hits, extract junctions, and toss the ones that arent strongly enough supported. |
1409 |
– |
|
1540 |
|
filter_junctions(junctions, gtf_junctions); |
1541 |
|
|
1542 |
|
//size_t num_juncs_after_filter = junctions.size(); |
1544 |
|
// num_unfiltered_juncs - num_juncs_after_filter); |
1545 |
|
|
1546 |
|
/* |
1547 |
< |
size_t small_overhangs = 0; |
1548 |
< |
for (JunctionSet::iterator i = junctions.begin(); i != junctions.end(); ++i) |
1547 |
> |
size_t small_overhangs = 0; |
1548 |
> |
for (JunctionSet::iterator i = junctions.begin(); i != junctions.end(); ++i) |
1549 |
|
{ |
1550 |
< |
if (i->second.accepted && |
1551 |
< |
(i->second.left_extent < min_anchor_len || i->second.right_extent < min_anchor_len)) |
1552 |
< |
{ |
1553 |
< |
small_overhangs++; |
1554 |
< |
} |
1550 |
> |
if (i->second.accepted && |
1551 |
> |
(i->second.left_extent < min_anchor_len || i->second.right_extent < min_anchor_len)) |
1552 |
> |
{ |
1553 |
> |
small_overhangs++; |
1554 |
> |
} |
1555 |
|
} |
1556 |
|
|
1557 |
< |
if (small_overhangs >0) |
1558 |
< |
fprintf(stderr, "Warning: %lu small overhang junctions!\n", (long unsigned int)small_overhangs); |
1557 |
> |
if (small_overhangs >0) |
1558 |
> |
fprintf(stderr, "Warning: %lu small overhang junctions!\n", (long unsigned int)small_overhangs); |
1559 |
|
*/ |
1560 |
|
|
1561 |
|
JunctionSet vfinal_junctions[num_threads]; |
1567 |
|
{ |
1568 |
|
ReportWorker worker; |
1569 |
|
|
1570 |
+ |
worker.sam_header_fname = sam_header; |
1571 |
|
char filename[1024] = {0}; |
1572 |
|
sprintf(filename, "%s%d.bam", bam_output_fname.c_str(), i); |
1573 |
|
worker.bam_output_fname = filename; |
1574 |
< |
worker.sam_header_fname = sam_header; |
1575 |
< |
sprintf(filename, "%s/unmapped_left%d.fq", output_dir.c_str(), i); |
1576 |
< |
worker.left_um_fname = filename; |
1574 |
> |
string tmpoutdir=getFdir(worker.bam_output_fname); |
1575 |
> |
worker.left_um_fname = tmpoutdir; |
1576 |
> |
sprintf(filename, "unmapped_left_%d.bam", i); |
1577 |
> |
worker.left_um_fname+=filename; |
1578 |
|
|
1579 |
|
if (right_reads_fname != "") |
1580 |
|
{ |
1581 |
< |
sprintf(filename, "%s/unmapped_right%d.fq", output_dir.c_str(), i); |
1582 |
< |
worker.right_um_fname = filename; |
1581 |
> |
sprintf(filename, "unmapped_right_%d.bam", i); |
1582 |
> |
worker.right_um_fname = tmpoutdir; |
1583 |
> |
worker.right_um_fname += filename; |
1584 |
|
} |
1585 |
|
|
1586 |
|
worker.left_reads_fname = left_reads_fname; |
1590 |
|
|
1591 |
|
worker.gtf_junctions = >f_junctions; |
1592 |
|
worker.junctions = &junctions; |
1593 |
+ |
worker.insertions = &insertions; |
1594 |
+ |
worker.deletions = &deletions; |
1595 |
|
worker.fusions = &fusions; |
1596 |
+ |
worker.coverage = &coverage; |
1597 |
+ |
|
1598 |
|
worker.final_junctions = &vfinal_junctions[i]; |
1599 |
|
worker.final_insertions = &vfinal_insertions[i]; |
1600 |
|
worker.final_deletions = &vfinal_deletions[i]; |