1 |
gpertea |
29 |
/* |
2 |
|
|
* bwt_map.cpp |
3 |
|
|
* TopHat |
4 |
|
|
* |
5 |
|
|
* Created by Cole Trapnell on 11/17/08. |
6 |
|
|
* Copyright 2008 Cole Trapnell. All rights reserved. |
7 |
|
|
* |
8 |
|
|
*/ |
9 |
|
|
|
10 |
|
|
#ifdef HAVE_CONFIG_H |
11 |
|
|
#include <config.h> |
12 |
|
|
#endif |
13 |
|
|
|
14 |
|
|
#include <cassert> |
15 |
|
|
#include <cstdlib> |
16 |
|
|
#include <iostream> |
17 |
|
|
#include <set> |
18 |
|
|
#include <vector> |
19 |
|
|
#include <cmath> |
20 |
|
|
|
21 |
|
|
#include "common.h" |
22 |
|
|
#include "bwt_map.h" |
23 |
|
|
#include "tokenize.h" |
24 |
|
|
#include "reads.h" |
25 |
|
|
|
26 |
|
|
using namespace std; |
27 |
|
|
|
28 |
|
|
void HitTable::add_hit(const BowtieHit& bh, bool check_uniqueness) |
29 |
|
|
{ |
30 |
|
|
uint32_t reference_id = bh.ref_id(); |
31 |
|
|
|
32 |
|
|
pair<RefHits::iterator, bool> ret = |
33 |
|
|
_hits_for_ref.insert(make_pair(reference_id, HitList())); |
34 |
|
|
HitList& hl = ret.first->second; |
35 |
|
|
|
36 |
|
|
if (check_uniqueness) |
37 |
|
|
{ |
38 |
|
|
// Check uniqueness, in case we are adding spliced hits from |
39 |
|
|
// several spliced alignment sources (e.g. de novo hashing + Bowtie |
40 |
|
|
// against a user-supplied index). We don't want to count the same |
41 |
|
|
// alignment twice if it happened to be found by more than one method |
42 |
|
|
HitList::const_iterator lb = lower_bound(hl.begin(), hl.end(), bh, hit_insert_id_lt); |
43 |
|
|
HitList::const_iterator ub = upper_bound(hl.begin(), hl.end(), bh, hit_insert_id_lt); |
44 |
|
|
for (; lb != ub && lb != hl.end(); ++lb) |
45 |
|
|
{ |
46 |
|
|
if (*lb == bh) |
47 |
|
|
{ |
48 |
|
|
//fprintf(stderr, "Chucking duplicate read %d by identity\n", bh.insert_id()); |
49 |
|
|
return; |
50 |
|
|
} |
51 |
|
|
|
52 |
|
|
if (lb->insert_id() == bh.insert_id() && |
53 |
|
|
lb->ref_id() == bh.ref_id() && |
54 |
|
|
lb->antisense_align() == bh.antisense_align()) |
55 |
|
|
{ |
56 |
|
|
// If we get here, we may be looking at the same alignment |
57 |
|
|
// However, spanning_reads may report a shorter, trimmed alignment |
58 |
|
|
// so not all fields will be equal. If they just disagree on the |
59 |
|
|
// ends, and don't indicate a different junction coord, the |
60 |
|
|
// alignments are the same. |
61 |
|
|
|
62 |
|
|
if ((lb->left() <= bh.left() && lb->right() >= bh.right()) || |
63 |
|
|
(bh.left() <= lb->left() && bh.right() >= lb->right())) |
64 |
|
|
{ |
65 |
|
|
vector<pair<int,int> > lb_gaps, bh_gaps; |
66 |
|
|
lb->gaps(lb_gaps); |
67 |
|
|
bh.gaps(bh_gaps); |
68 |
|
|
if (lb_gaps == bh_gaps) |
69 |
|
|
{ |
70 |
|
|
// One alignment is contained in the other, they agree on |
71 |
|
|
// where the gaps, if any, are, and they share an id |
72 |
|
|
// => this is a redundant aligment, so toss it |
73 |
|
|
//fprintf(stderr, "Chucking duplicate read %d by gap agreement\n", bh.insert_id()); |
74 |
|
|
return; |
75 |
|
|
} |
76 |
|
|
} |
77 |
|
|
} |
78 |
|
|
} |
79 |
|
|
} |
80 |
|
|
_total_hits++; |
81 |
|
|
hl.push_back(bh); |
82 |
|
|
} |
83 |
|
|
|
84 |
|
|
bool hit_insert_id_lt(const BowtieHit& h1, const BowtieHit& h2) |
85 |
|
|
{ |
86 |
|
|
return h1.insert_id() < h2.insert_id(); |
87 |
|
|
} |
88 |
|
|
|
89 |
|
|
void LineHitFactory::openStream(HitStream& hs) |
90 |
|
|
{ |
91 |
|
|
if (hs._hit_file==NULL && !hs._hit_file_name.empty()) { |
92 |
|
|
//open the file for HitStream here |
93 |
|
|
hs._hit_file=fopen(hs._hit_file_name.c_str(),"r"); |
94 |
|
|
if (hs._hit_file==NULL) |
95 |
|
|
err_die("Error opening HitStream file %s\n",hs._hit_file_name.c_str()); |
96 |
|
|
return; |
97 |
|
|
} |
98 |
|
|
if (hs._fzpipe!=NULL) { |
99 |
|
|
hs._hit_file=hs._fzpipe->file; |
100 |
|
|
} |
101 |
|
|
} |
102 |
|
|
|
103 |
|
|
void LineHitFactory::rewind(HitStream& hs) |
104 |
|
|
{ |
105 |
|
|
if (hs._fzpipe!=NULL) { |
106 |
|
|
hs._fzpipe->rewind(); |
107 |
|
|
hs._hit_file=hs._fzpipe->file; |
108 |
|
|
} |
109 |
|
|
else if (hs._hit_file) ::rewind((FILE*)(hs._hit_file)); |
110 |
|
|
} |
111 |
|
|
|
112 |
|
|
bool LineHitFactory::next_record(HitStream& hs, const char*& buf, size_t& buf_size) { |
113 |
|
|
FILE* f=(FILE *)(hs._hit_file); |
114 |
|
|
bool new_rec = (fgets(_hit_buf, _hit_buf_max_sz - 1, f)!=NULL); |
115 |
|
|
if (!new_rec || feof(f)) { |
116 |
|
|
hs._eof=true; |
117 |
|
|
return false; |
118 |
|
|
} |
119 |
|
|
++_line_num; |
120 |
|
|
char* nl = strrchr(_hit_buf, '\n'); |
121 |
|
|
if (nl) *nl = 0; |
122 |
|
|
buf = _hit_buf; |
123 |
|
|
buf_size = _hit_buf_max_sz - 1; |
124 |
|
|
return true; |
125 |
|
|
} |
126 |
|
|
|
127 |
|
|
void LineHitFactory::closeStream(HitStream& hs) { |
128 |
|
|
if (hs._fzpipe!=NULL) { |
129 |
|
|
hs._fzpipe->close(); |
130 |
|
|
return; |
131 |
|
|
} |
132 |
|
|
if (hs._hit_file!=NULL) { |
133 |
|
|
fclose((FILE*)(hs._hit_file)); |
134 |
|
|
hs._hit_file=NULL; |
135 |
|
|
} |
136 |
|
|
} |
137 |
|
|
void BAMHitFactory::openStream(HitStream& hs) { |
138 |
|
|
if (hs._hit_file==NULL) { |
139 |
|
|
if (hs._hit_file_name.empty()) |
140 |
|
|
//err_die("Error: invalid HitStream set for BAMHitFactory(file name missing)\n"); |
141 |
|
|
return; //invalid stream, could be just a place holder |
142 |
|
|
//open the file here if not already open |
143 |
|
|
string fext=getFext(hs._hit_file_name); |
144 |
|
|
if (fext=="sam") |
145 |
|
|
hs._hit_file = samopen(hs._hit_file_name.c_str(), "r", 0); |
146 |
|
|
else |
147 |
|
|
hs._hit_file = samopen(hs._hit_file_name.c_str(), "rb", 0); |
148 |
|
|
|
149 |
|
|
samfile_t* sam_file=(samfile_t*)(hs._hit_file); |
150 |
|
|
|
151 |
|
|
if (sam_file == NULL) |
152 |
|
|
err_die("Error opening SAM file %s\n", hs._hit_file_name.c_str()); |
153 |
|
|
if (sam_file->header == NULL) |
154 |
|
|
err_die("Error: no SAM header found for file %s\n", hs._hit_file_name.c_str()); |
155 |
|
|
memset(&_next_hit, 0, sizeof(_next_hit)); |
156 |
|
|
//_beginning = bgzf_tell(sam_file->x.bam); |
157 |
|
|
_sam_header=sam_file->header; |
158 |
|
|
if (inspect_header(hs) == false) |
159 |
|
|
err_die("Error: invalid SAM header for file %s\n", |
160 |
|
|
hs._hit_file_name.c_str()); |
161 |
|
|
} |
162 |
|
|
} |
163 |
|
|
|
164 |
|
|
void BAMHitFactory::closeStream(HitStream& hs) { |
165 |
|
|
if (hs._hit_file) { |
166 |
|
|
samclose((samfile_t*)(hs._hit_file)); |
167 |
|
|
} |
168 |
|
|
hs._hit_file=NULL; |
169 |
|
|
_sam_header=NULL; |
170 |
|
|
} |
171 |
|
|
|
172 |
|
|
void BAMHitFactory::rewind(HitStream& hs) |
173 |
|
|
{ |
174 |
|
|
/* |
175 |
|
|
if (_hit_file && ((samfile_t*)_hit_file)->x.bam) |
176 |
|
|
{ |
177 |
|
|
bgzf_seek(((samfile_t*)_hit_file)->x.bam, _beginning, SEEK_SET); |
178 |
|
|
_eof = false; |
179 |
|
|
} |
180 |
|
|
*/ |
181 |
|
|
this->closeStream(hs); |
182 |
|
|
this->openStream(hs); |
183 |
|
|
} |
184 |
|
|
|
185 |
|
|
string BAMHitFactory::hitfile_rec(HitStream& hs, const char* hit_buf) { |
186 |
|
|
const bam1_t* bamrec=(const bam1_t*)hit_buf; |
187 |
|
|
char* tamline=bam_format1(((samfile_t*)(hs._hit_file))->header, bamrec); |
188 |
|
|
string sam_line(tamline); |
189 |
|
|
free(tamline); |
190 |
|
|
return sam_line; |
191 |
|
|
} |
192 |
|
|
|
193 |
|
|
bool BAMHitFactory::next_record(HitStream& hs, const char*& buf, size_t& buf_size) { |
194 |
|
|
if (_next_hit.data) { |
195 |
|
|
free(_next_hit.data); |
196 |
|
|
_next_hit.data = NULL; |
197 |
|
|
} |
198 |
|
|
_sam_header=((samfile_t*)(hs._hit_file))->header; //needed by get_hit_from_buf later on |
199 |
|
|
if (hs.eof() || !hs.ready()) return false; |
200 |
|
|
|
201 |
|
|
//mark_curr_pos(); |
202 |
|
|
|
203 |
|
|
memset(&_next_hit, 0, sizeof(_next_hit)); |
204 |
|
|
|
205 |
|
|
int bytes_read = samread((samfile_t*)(hs._hit_file), &_next_hit); |
206 |
|
|
if (bytes_read <= 0) { |
207 |
|
|
hs._eof = true; |
208 |
|
|
return false; |
209 |
|
|
} |
210 |
|
|
buf = (const char*)&_next_hit; |
211 |
|
|
buf_size = bytes_read; |
212 |
|
|
return true; |
213 |
|
|
} |
214 |
|
|
|
215 |
|
|
|
216 |
|
|
BowtieHit HitFactory::create_hit(const string& insert_name, |
217 |
|
|
const string& ref_name, |
218 |
|
|
int left, |
219 |
|
|
const vector<CigarOp>& cigar, |
220 |
|
|
bool antisense_aln, |
221 |
|
|
bool antisense_splice, |
222 |
|
|
unsigned char edit_dist, |
223 |
|
|
unsigned char splice_mms, |
224 |
|
|
bool end) |
225 |
|
|
{ |
226 |
|
|
uint64_t insert_id = _insert_table.get_id(insert_name); |
227 |
|
|
uint32_t reference_id = _ref_table.get_id(ref_name, NULL, 0); |
228 |
|
|
|
229 |
|
|
return BowtieHit(reference_id, |
230 |
|
|
insert_id, |
231 |
|
|
left, |
232 |
|
|
cigar, |
233 |
|
|
antisense_aln, |
234 |
|
|
antisense_splice, |
235 |
|
|
edit_dist, |
236 |
|
|
splice_mms, |
237 |
|
|
end); |
238 |
|
|
} |
239 |
|
|
|
240 |
|
|
BowtieHit HitFactory::create_hit(const string& insert_name, |
241 |
|
|
const string& ref_name, |
242 |
|
|
uint32_t left, |
243 |
|
|
uint32_t read_len, |
244 |
|
|
bool antisense_aln, |
245 |
|
|
unsigned char edit_dist, |
246 |
|
|
bool end) |
247 |
|
|
{ |
248 |
|
|
uint64_t insert_id = _insert_table.get_id(insert_name); |
249 |
|
|
uint32_t reference_id = _ref_table.get_id(ref_name, NULL, 0); |
250 |
|
|
|
251 |
|
|
return BowtieHit(reference_id, |
252 |
|
|
insert_id, |
253 |
|
|
left, |
254 |
|
|
read_len, |
255 |
|
|
antisense_aln, |
256 |
|
|
edit_dist, |
257 |
|
|
end); |
258 |
|
|
} |
259 |
|
|
|
260 |
|
|
bool BowtieHitFactory::get_hit_from_buf(const char* orig_bwt_buf, |
261 |
|
|
BowtieHit& bh, |
262 |
|
|
bool strip_slash, |
263 |
|
|
char* name_out, |
264 |
|
|
char* name_tags, |
265 |
|
|
char* seq, |
266 |
|
|
char* qual) |
267 |
|
|
{ |
268 |
|
|
if (!orig_bwt_buf || !*orig_bwt_buf) |
269 |
|
|
return false; |
270 |
|
|
|
271 |
|
|
static const int buf_size = 2048; |
272 |
|
|
|
273 |
|
|
char bwt_buf[buf_size]; |
274 |
|
|
strcpy(bwt_buf, orig_bwt_buf); |
275 |
|
|
//const char* bwt_fmt_str = "%s %c %s %d %s %s %d %s"; |
276 |
|
|
|
277 |
|
|
char orientation; |
278 |
|
|
|
279 |
|
|
//int bwtf_ret = 0; |
280 |
|
|
//uint32_t seqid = 0; |
281 |
|
|
|
282 |
|
|
char text_name[buf_size]; |
283 |
|
|
unsigned int text_offset; |
284 |
|
|
|
285 |
|
|
|
286 |
|
|
//unsigned int other_occs; |
287 |
|
|
char mismatches[buf_size]; |
288 |
|
|
//memset(mismatches, 0, sizeof(mismatches)); |
289 |
|
|
|
290 |
|
|
const char* buf = bwt_buf; |
291 |
|
|
char* name = get_token((char**)&buf,"\t"); |
292 |
|
|
char* orientation_str = get_token((char**)&buf,"\t"); |
293 |
|
|
char* text_name_str = get_token((char**)&buf,"\t"); |
294 |
|
|
|
295 |
|
|
strcpy(text_name, text_name_str); |
296 |
|
|
|
297 |
|
|
char* text_offset_str = get_token((char**)&buf,"\t"); |
298 |
|
|
char* seq_str = get_token((char**)&buf,"\t"); |
299 |
|
|
if (seq) |
300 |
|
|
strcpy(seq, seq_str); |
301 |
|
|
|
302 |
|
|
const char* qual_str = get_token((char**)&buf,"\t"); |
303 |
|
|
if (qual) |
304 |
|
|
strcpy(qual, qual_str); |
305 |
|
|
|
306 |
|
|
/*const char* other_occs_str =*/ get_token((char**)&buf, "\t"); |
307 |
|
|
mismatches[0] = 0; |
308 |
|
|
char* mismatches_str = get_token((char**)&buf, "\t"); |
309 |
|
|
if (mismatches_str) |
310 |
|
|
strcpy(mismatches, mismatches_str); |
311 |
|
|
|
312 |
|
|
orientation = orientation_str[0]; |
313 |
|
|
text_offset = atoi(text_offset_str); |
314 |
|
|
|
315 |
|
|
bool end = true; |
316 |
|
|
unsigned int seg_offset = 0; |
317 |
|
|
unsigned int seg_num = 0; |
318 |
|
|
unsigned int num_segs = 0; |
319 |
|
|
|
320 |
|
|
// Copy the tag out of the name field before we might wipe it out |
321 |
|
|
char* pipe = strrchr(name, '|'); |
322 |
|
|
if (pipe) |
323 |
|
|
{ |
324 |
|
|
if (name_tags) |
325 |
|
|
strcpy(name_tags, pipe); |
326 |
|
|
|
327 |
|
|
char* tag_buf = pipe + 1; |
328 |
|
|
if (strchr(tag_buf, ':')) |
329 |
|
|
{ |
330 |
|
|
sscanf(tag_buf, "%u:%u:%u", &seg_offset, &seg_num, &num_segs); |
331 |
|
|
if (seg_num + 1 == num_segs) |
332 |
|
|
end = true; |
333 |
|
|
else |
334 |
|
|
end = false; |
335 |
|
|
} |
336 |
|
|
|
337 |
|
|
*pipe = 0; |
338 |
|
|
} |
339 |
|
|
// Stripping the slash and number following it gives the insert name |
340 |
|
|
char* slash = strrchr(name, '/'); |
341 |
|
|
if (strip_slash && slash) |
342 |
|
|
*slash = 0; |
343 |
|
|
|
344 |
|
|
size_t read_len = strlen(seq_str); |
345 |
|
|
|
346 |
|
|
// Add this alignment to the table of hits for this half of the |
347 |
|
|
// Bowtie map |
348 |
|
|
|
349 |
|
|
//vector<string> mismatch_toks; |
350 |
|
|
char* pch = strtok (mismatches,","); |
351 |
|
|
unsigned char num_mismatches = 0; |
352 |
|
|
while (pch != NULL) |
353 |
|
|
{ |
354 |
|
|
char* colon = strchr(pch, ':'); |
355 |
|
|
if (colon) |
356 |
|
|
{ |
357 |
|
|
num_mismatches++; |
358 |
|
|
} |
359 |
|
|
//mismatch_toks.push_back(pch); |
360 |
|
|
pch = strtok (NULL, ","); |
361 |
|
|
} |
362 |
|
|
|
363 |
|
|
bh = create_hit(name, |
364 |
|
|
text_name, |
365 |
|
|
text_offset, |
366 |
|
|
(int)read_len, |
367 |
|
|
orientation == '-', |
368 |
|
|
num_mismatches, |
369 |
|
|
end); |
370 |
|
|
|
371 |
|
|
return true; |
372 |
|
|
} |
373 |
|
|
|
374 |
|
|
int anchor_mismatch = 0; |
375 |
|
|
|
376 |
gpertea |
70 |
|
377 |
|
|
void parseSegReadName(char* name, char*& name_tags, bool strip_slash, |
378 |
|
|
bool &end, unsigned int &seg_offset, unsigned int& seg_num, |
379 |
|
|
unsigned int & num_segs) { |
380 |
|
|
char* pipe = strrchr(name, '|'); |
381 |
|
|
if (pipe) |
382 |
|
|
{ |
383 |
|
|
if (name_tags) |
384 |
|
|
strcpy(name_tags, pipe); |
385 |
|
|
|
386 |
|
|
char* tag_buf = pipe + 1; |
387 |
|
|
if (strchr(tag_buf, ':')) |
388 |
|
|
{ |
389 |
|
|
sscanf(tag_buf, "%u:%u:%u", &seg_offset, &seg_num, &num_segs); |
390 |
|
|
if (seg_num + 1 == num_segs) |
391 |
|
|
end = true; |
392 |
|
|
else |
393 |
|
|
end = false; |
394 |
|
|
} |
395 |
|
|
|
396 |
|
|
*pipe = 0; |
397 |
|
|
} |
398 |
|
|
|
399 |
|
|
// Stripping the slash and number following it gives the insert name |
400 |
|
|
char* slash = strrchr(name, '/'); |
401 |
|
|
if (strip_slash && slash) |
402 |
|
|
*slash = 0; |
403 |
|
|
} |
404 |
|
|
|
405 |
|
|
|
406 |
gpertea |
29 |
bool SplicedBowtieHitFactory::get_hit_from_buf(const char* orig_bwt_buf, |
407 |
|
|
BowtieHit& bh, |
408 |
|
|
bool strip_slash, |
409 |
|
|
char* name_out, |
410 |
|
|
char* name_tags, |
411 |
|
|
char* seq, |
412 |
|
|
char* qual) |
413 |
|
|
{ |
414 |
|
|
if (!orig_bwt_buf || !*orig_bwt_buf) |
415 |
|
|
return false; |
416 |
|
|
|
417 |
|
|
static const int buf_size = 2048; |
418 |
|
|
|
419 |
|
|
char bwt_buf[buf_size]; |
420 |
|
|
strcpy(bwt_buf, orig_bwt_buf); |
421 |
|
|
//const char* bwt_fmt_str = "%s %c %s %d %s %s %d %s"; |
422 |
|
|
|
423 |
|
|
char orientation; |
424 |
|
|
char text_name[buf_size]; |
425 |
|
|
unsigned int text_offset; |
426 |
|
|
char mismatches[buf_size]; |
427 |
|
|
//memset(mismatches, 0, sizeof(mismatches)); |
428 |
|
|
|
429 |
|
|
char* buf = bwt_buf; |
430 |
|
|
char* name = get_token((char**)&buf,"\t"); |
431 |
|
|
char* orientation_str = get_token((char**)&buf,"\t"); |
432 |
|
|
char* text_name_str = get_token((char**)&buf,"\t"); |
433 |
|
|
strcpy(text_name, text_name_str); |
434 |
|
|
|
435 |
|
|
char* text_offset_str = get_token((char**)&buf,"\t"); |
436 |
|
|
char* seq_str = get_token((char**)&buf,"\t"); |
437 |
|
|
if (seq) |
438 |
|
|
strcpy(seq, seq_str); |
439 |
|
|
|
440 |
|
|
const char* qual_str = get_token((char**)&buf,"\t"); |
441 |
|
|
if (qual) |
442 |
|
|
strcpy(qual, qual_str); |
443 |
|
|
|
444 |
|
|
/*const char* other_occs_str =*/ get_token((char**)&buf, "\t"); |
445 |
|
|
mismatches[0] = 0; |
446 |
|
|
char* mismatches_str = get_token((char**)&buf, "\t"); |
447 |
|
|
if (mismatches_str) |
448 |
|
|
strcpy(mismatches, mismatches_str); |
449 |
|
|
|
450 |
|
|
orientation = orientation_str[0]; |
451 |
|
|
text_offset = atoi(text_offset_str); |
452 |
|
|
|
453 |
|
|
bool end = true; |
454 |
|
|
unsigned int seg_offset = 0; |
455 |
|
|
unsigned int seg_num = 0; |
456 |
|
|
unsigned int num_segs = 0; |
457 |
|
|
|
458 |
|
|
// Copy the tag out of the name field before we might wipe it out |
459 |
gpertea |
70 |
parseSegReadName(name, name_tags, strip_slash, end, seg_offset, seg_num, num_segs); |
460 |
gpertea |
29 |
|
461 |
|
|
// Add this alignment to the table of hits for this half of the |
462 |
|
|
// Bowtie map |
463 |
|
|
|
464 |
|
|
// Parse the text_name field to recover the splice coords |
465 |
|
|
vector<string> toks; |
466 |
|
|
|
467 |
|
|
tokenize_strict(text_name, "|", toks); |
468 |
|
|
|
469 |
|
|
int num_extra_toks = (int)toks.size() - 6; |
470 |
|
|
|
471 |
|
|
if (num_extra_toks >= 0) |
472 |
|
|
{ |
473 |
|
|
static const uint8_t left_window_edge_field = 1; |
474 |
|
|
static const uint8_t splice_field = 2; |
475 |
|
|
//static const uint8_t right_window_edge_field = 3; |
476 |
|
|
static const uint8_t junction_type_field = 4; |
477 |
|
|
static const uint8_t strand_field = 5; |
478 |
|
|
|
479 |
|
|
string contig = toks[0]; |
480 |
|
|
for (int t = 1; t <= num_extra_toks; ++t) |
481 |
|
|
{ |
482 |
|
|
contig += "|"; |
483 |
|
|
contig += toks[t]; |
484 |
|
|
} |
485 |
|
|
|
486 |
|
|
vector<string> splice_toks; |
487 |
|
|
tokenize(toks[num_extra_toks + splice_field], "-", splice_toks); |
488 |
|
|
if (splice_toks.size() != 2) |
489 |
|
|
{ |
490 |
|
|
fprintf(stderr, "Warning: found malformed splice record, skipping:\n"); |
491 |
|
|
//fprintf(stderr, "%s (token: %s)\n", text_name, |
492 |
|
|
// toks[num_extra_toks + splice_field].c_str()); |
493 |
|
|
return false; |
494 |
|
|
} |
495 |
|
|
|
496 |
|
|
// |
497 |
|
|
// check for an insertion hit |
498 |
|
|
// |
499 |
|
|
if(toks[num_extra_toks + junction_type_field] == "ins") |
500 |
|
|
{ |
501 |
|
|
int8_t spliced_read_len = strlen(seq_str); |
502 |
|
|
/* |
503 |
|
|
* The 0-based position of the left edge of the alignment. Note that this |
504 |
|
|
* value may need to be futher corrected to account for the presence of |
505 |
|
|
* of the insertion. |
506 |
|
|
*/ |
507 |
|
|
uint32_t left = atoi(toks[num_extra_toks + left_window_edge_field].c_str()) + text_offset; |
508 |
|
|
uint32_t right = left + spliced_read_len - 1; |
509 |
|
|
|
510 |
gpertea |
137 |
|
511 |
gpertea |
29 |
/* |
512 |
|
|
* The 0-based position of the last genomic sequence before the insertion |
513 |
|
|
*/ |
514 |
|
|
uint32_t left_splice_pos = atoi(splice_toks[0].c_str()); |
515 |
|
|
|
516 |
|
|
string insertedSequence = splice_toks[1]; |
517 |
|
|
/* |
518 |
|
|
* The 0-based position of the first genomic sequence after teh insertion |
519 |
|
|
*/ |
520 |
|
|
uint32_t right_splice_pos = left_splice_pos + 1; |
521 |
|
|
if(left > left_splice_pos){ |
522 |
|
|
/* |
523 |
|
|
* The genomic position of the left edge of the alignment needs to be corrected |
524 |
|
|
* If the alignment does not extend into the insertion, simply subtract the length |
525 |
|
|
* of the inserted sequence, otherwise, just set it equal to the right edge |
526 |
|
|
*/ |
527 |
|
|
left = left - insertedSequence.length(); |
528 |
|
|
if(left < right_splice_pos){ |
529 |
|
|
left = right_splice_pos; |
530 |
|
|
} |
531 |
|
|
} |
532 |
|
|
if(right > left_splice_pos){ |
533 |
|
|
right = right - insertedSequence.length(); |
534 |
|
|
if(right < left_splice_pos){ |
535 |
|
|
right = left_splice_pos; |
536 |
|
|
} |
537 |
|
|
} |
538 |
|
|
/* |
539 |
|
|
* Now, right and left should be properly transformed into genomic coordinates |
540 |
|
|
* We should be able to deduce how much the alignment matches the insertion |
541 |
|
|
* simply based on the length of the read |
542 |
|
|
*/ |
543 |
|
|
int left_match_length = 0; |
544 |
|
|
if(left <= left_splice_pos){ |
545 |
|
|
left_match_length = left_splice_pos - left + 1; |
546 |
|
|
} |
547 |
|
|
int right_match_length = 0; |
548 |
|
|
if(right >= right_splice_pos){ |
549 |
|
|
right_match_length = right - right_splice_pos + 1; |
550 |
|
|
} |
551 |
|
|
int insertion_match_length = spliced_read_len - left_match_length - right_match_length; |
552 |
|
|
|
553 |
|
|
if(left_match_length <= 0 || right_match_length <= 0 || insertion_match_length <= 0) |
554 |
|
|
return false; |
555 |
|
|
|
556 |
|
|
string junction_strand = toks[num_extra_toks + strand_field]; |
557 |
|
|
if(junction_strand != "rev" && junction_strand != "fwd"){ |
558 |
|
|
fprintf(stderr, "Malformed insertion record\n"); |
559 |
|
|
return false; |
560 |
|
|
} |
561 |
|
|
|
562 |
|
|
char* pch = strtok( mismatches, ","); |
563 |
|
|
unsigned char num_mismatches = 0; |
564 |
|
|
|
565 |
|
|
/* |
566 |
|
|
* remember that text_offset holds the left end of the |
567 |
|
|
* alignment, relative to the start of the contig |
568 |
|
|
*/ |
569 |
|
|
|
570 |
|
|
/* |
571 |
|
|
* The 0-based relative position of the left-most character |
572 |
|
|
* before the insertion in the contig |
573 |
|
|
*/ |
574 |
|
|
int relative_splice_pos = left_splice_pos - atoi(toks[num_extra_toks + left_window_edge_field].c_str()); |
575 |
|
|
while (pch != NULL) |
576 |
|
|
{ |
577 |
|
|
char* colon = strchr(pch, ':'); |
578 |
|
|
if (colon) |
579 |
|
|
{ |
580 |
|
|
*colon = 0; |
581 |
|
|
int mismatch_pos = atoi(pch); |
582 |
|
|
|
583 |
|
|
/* |
584 |
|
|
* for reversely mapped reads, |
585 |
|
|
* find the correct mismatched position. |
586 |
|
|
*/ |
587 |
|
|
if(orientation == '-'){ |
588 |
|
|
mismatch_pos = spliced_read_len - mismatch_pos - 1; |
589 |
|
|
} |
590 |
|
|
|
591 |
|
|
/* |
592 |
|
|
* Only count mismatches outside of the insertion region |
593 |
|
|
* If there is a mismatch within the insertion, |
594 |
|
|
* disallow this hit |
595 |
|
|
*/ |
596 |
|
|
if(mismatch_pos + text_offset <= relative_splice_pos || mismatch_pos + text_offset > relative_splice_pos + insertedSequence.length()){ |
597 |
|
|
num_mismatches++; |
598 |
|
|
}else{ |
599 |
|
|
return false; |
600 |
|
|
} |
601 |
|
|
} |
602 |
|
|
pch = strtok (NULL, ","); |
603 |
|
|
} |
604 |
|
|
|
605 |
|
|
|
606 |
|
|
vector<CigarOp> cigar; |
607 |
|
|
cigar.push_back(CigarOp(MATCH, left_match_length)); |
608 |
|
|
cigar.push_back(CigarOp(INS, insertion_match_length)); |
609 |
|
|
cigar.push_back(CigarOp(MATCH, right_match_length)); |
610 |
|
|
|
611 |
|
|
/* |
612 |
|
|
* For now, disallow hits that don't span |
613 |
|
|
* the insertion. If we allow these types of hits, |
614 |
|
|
* then long_spanning.cpp needs to be updated |
615 |
|
|
* in order to intelligently merge these kinds |
616 |
|
|
* of reads back together |
617 |
|
|
* |
618 |
|
|
* Following code has been changed to allow segment that end |
619 |
|
|
* in an insertion |
620 |
|
|
*/ |
621 |
|
|
bh = create_hit(name, |
622 |
|
|
contig, |
623 |
|
|
left, |
624 |
|
|
cigar, |
625 |
|
|
orientation == '-', |
626 |
|
|
junction_strand == "rev", |
627 |
|
|
num_mismatches, |
628 |
|
|
0, |
629 |
|
|
end); |
630 |
|
|
return true; |
631 |
|
|
} |
632 |
|
|
|
633 |
|
|
else |
634 |
|
|
{ |
635 |
|
|
uint32_t left = atoi(toks[num_extra_toks + left_window_edge_field].c_str()) + text_offset; |
636 |
|
|
int spliced_read_len = strlen(seq_str); |
637 |
|
|
int8_t left_splice_pos = atoi(splice_toks[0].c_str()) - left + 1; |
638 |
|
|
if(left_splice_pos > spliced_read_len) left_splice_pos = spliced_read_len; |
639 |
|
|
int8_t right_splice_pos = spliced_read_len - left_splice_pos; |
640 |
|
|
|
641 |
|
|
if (right_splice_pos <= 0 || left_splice_pos <= 0) |
642 |
|
|
return false; |
643 |
|
|
|
644 |
|
|
if (orientation == '+') |
645 |
|
|
{ |
646 |
|
|
if (left_splice_pos + seg_offset < _anchor_length){ |
647 |
|
|
return false; |
648 |
|
|
} |
649 |
|
|
} |
650 |
|
|
else |
651 |
|
|
{ |
652 |
|
|
if (right_splice_pos + seg_offset < _anchor_length) |
653 |
|
|
return false; |
654 |
|
|
} |
655 |
|
|
//uint32_t right = atoi(splice_toks[1].c_str()) + right_splice_pos; |
656 |
|
|
//atoi(toks[right_window_edge_field].c_str()); |
657 |
|
|
int gap_len = atoi(splice_toks[1].c_str()) - atoi(splice_toks[0].c_str()) - 1; |
658 |
|
|
|
659 |
|
|
string junction_strand = toks[num_extra_toks + strand_field]; |
660 |
|
|
if (!(junction_strand == "rev" || junction_strand == "fwd")|| |
661 |
|
|
!(orientation == '-' || orientation == '+')) |
662 |
|
|
{ |
663 |
|
|
fprintf(stderr, "Warning: found malformed splice record, skipping\n"); |
664 |
|
|
//fprintf(stderr, "junction_strand=%s, orientation='%c'\n", |
665 |
|
|
// junction_strand.c_str(), orientation); |
666 |
|
|
return false; |
667 |
|
|
} |
668 |
|
|
|
669 |
|
|
//vector<string> mismatch_toks; |
670 |
|
|
char* pch = strtok (mismatches,","); |
671 |
|
|
int mismatches_in_anchor = 0; |
672 |
|
|
unsigned char num_mismatches = 0; |
673 |
|
|
while (pch != NULL) |
674 |
|
|
{ |
675 |
|
|
char* colon = strchr(pch, ':'); |
676 |
|
|
if (colon) |
677 |
|
|
{ |
678 |
|
|
*colon = 0; |
679 |
|
|
num_mismatches++; |
680 |
|
|
int mismatch_pos = atoi(pch); |
681 |
|
|
if ((orientation == '+' && abs(mismatch_pos - left_splice_pos) < (int)min_anchor_len) || |
682 |
|
|
(orientation == '-' && abs(((int)spliced_read_len - left_splice_pos + 1) - mismatch_pos)) < (int)min_anchor_len) |
683 |
|
|
mismatches_in_anchor++; |
684 |
|
|
} |
685 |
|
|
//mismatch_toks.push_back(pch); |
686 |
|
|
pch = strtok (NULL, ","); |
687 |
|
|
} |
688 |
|
|
|
689 |
|
|
// FIXME: we probably should exclude these hits somewhere, but this |
690 |
|
|
// isn't the right place |
691 |
|
|
vector<CigarOp> cigar; |
692 |
|
|
cigar.push_back(CigarOp(MATCH, left_splice_pos)); |
693 |
|
|
if(toks[num_extra_toks + junction_type_field] == "del"){ |
694 |
|
|
cigar.push_back(CigarOp(DEL, gap_len)); |
695 |
|
|
}else{ |
696 |
|
|
cigar.push_back(CigarOp(REF_SKIP, gap_len)); |
697 |
|
|
} |
698 |
|
|
cigar.push_back(CigarOp(MATCH, right_splice_pos)); |
699 |
|
|
|
700 |
|
|
bh = create_hit(name, |
701 |
|
|
contig, |
702 |
|
|
left, |
703 |
|
|
cigar, |
704 |
|
|
orientation == '-', |
705 |
|
|
junction_strand == "rev", |
706 |
|
|
num_mismatches, |
707 |
|
|
mismatches_in_anchor, |
708 |
|
|
end); |
709 |
|
|
return true; |
710 |
|
|
} |
711 |
|
|
} |
712 |
|
|
else |
713 |
|
|
{ |
714 |
|
|
fprintf(stderr, "Warning: found malformed splice record, skipping\n"); |
715 |
|
|
//fprintf(stderr, "%s\n", orig_bwt_buf); |
716 |
|
|
// continue; |
717 |
|
|
return false; |
718 |
|
|
} |
719 |
|
|
|
720 |
|
|
return false; |
721 |
|
|
} |
722 |
|
|
|
723 |
gpertea |
70 |
|
724 |
gpertea |
137 |
int parseCigar(vector<CigarOp>& cigar, const char* cigar_str, |
725 |
gpertea |
70 |
bool &spliced_alignment) { |
726 |
|
|
const char* p_cig = cigar_str; |
727 |
gpertea |
137 |
int refspan=0; //alignment span on reference sequence |
728 |
gpertea |
70 |
while (*p_cig) |
729 |
|
|
{ |
730 |
|
|
char* t; |
731 |
gpertea |
137 |
int op_len = (int)strtol(p_cig, &t, 10); |
732 |
|
|
if (op_len <= 0) |
733 |
gpertea |
70 |
{ |
734 |
gpertea |
137 |
fprintf (stderr, "Error: CIGAR op has zero length\n"); |
735 |
|
|
return 0; |
736 |
gpertea |
70 |
} |
737 |
|
|
char op_char = toupper(*t); |
738 |
|
|
CigarOpCode opcode; |
739 |
gpertea |
137 |
switch (op_char) { |
740 |
|
|
case '=': |
741 |
|
|
case 'X': |
742 |
|
|
case 'M': opcode = MATCH; |
743 |
|
|
refspan+=op_len; |
744 |
|
|
break; |
745 |
|
|
case 'I': opcode = INS; |
746 |
|
|
break; |
747 |
|
|
case 'D': opcode = DEL; |
748 |
|
|
refspan+=op_len; |
749 |
|
|
break; |
750 |
|
|
case 'N': if (op_len > max_report_intron_length) |
751 |
|
|
return 0; |
752 |
|
|
opcode = REF_SKIP; |
753 |
|
|
spliced_alignment = true; |
754 |
|
|
refspan+=op_len; |
755 |
|
|
break; |
756 |
|
|
case 'S': opcode = SOFT_CLIP; |
757 |
|
|
break; |
758 |
|
|
case 'H': opcode = HARD_CLIP; |
759 |
|
|
break; |
760 |
|
|
case 'P': opcode = PAD; |
761 |
|
|
break; |
762 |
|
|
default: fprintf (stderr, "Error: invalid CIGAR operation\n"); |
763 |
|
|
return 0; |
764 |
|
|
} |
765 |
gpertea |
70 |
p_cig = t + 1; |
766 |
gpertea |
137 |
cigar.push_back(CigarOp(opcode, op_len)); |
767 |
|
|
} //while cigar codes |
768 |
gpertea |
70 |
if (*p_cig) { |
769 |
gpertea |
137 |
fprintf (stderr, "Error: unmatched CIGAR operation (%s)\n",p_cig); |
770 |
|
|
return 0; |
771 |
gpertea |
70 |
} |
772 |
gpertea |
137 |
return refspan; |
773 |
gpertea |
70 |
} |
774 |
|
|
|
775 |
gpertea |
137 |
int getSAMmismatches(char* &buf, vector<CigarOp>& cigar, |
776 |
|
|
vector<int>& mismatches, int& sam_nm, bool& antisense_splice) { |
777 |
|
|
int gspan=0;//genomic span of the alignment |
778 |
|
|
const char* tag_buf = buf; |
779 |
gpertea |
70 |
sam_nm=0; |
780 |
gpertea |
137 |
//num_mismatches=0; |
781 |
gpertea |
70 |
while((tag_buf = get_token((char**)&buf,"\t"))) |
782 |
|
|
{ |
783 |
|
|
vector<string> tuple_fields; |
784 |
|
|
tokenize(tag_buf,":", tuple_fields); |
785 |
|
|
if (tuple_fields.size() == 3) |
786 |
|
|
{ |
787 |
|
|
if (tuple_fields[0] == "XS") |
788 |
|
|
{ |
789 |
|
|
if (tuple_fields[2] == "-") |
790 |
|
|
antisense_splice = true; |
791 |
|
|
} |
792 |
|
|
else if (tuple_fields[0] == "NM") |
793 |
|
|
{ |
794 |
|
|
sam_nm = atoi(tuple_fields[2].c_str()); |
795 |
|
|
} |
796 |
|
|
else if (tuple_fields[0] == "NS") |
797 |
|
|
{ |
798 |
|
|
//ignored for now |
799 |
|
|
} |
800 |
|
|
else if (tuple_fields[0] == "MD") |
801 |
|
|
{ |
802 |
|
|
const char* p=tuple_fields[2].c_str(); |
803 |
|
|
int bi=0; //base offset position in the read |
804 |
|
|
while (*p != 0) { |
805 |
|
|
if (isdigit(*p)) { |
806 |
|
|
int v=atoi(p); |
807 |
|
|
do { p++; } while (isdigit(*p)); |
808 |
|
|
bi+=v; |
809 |
|
|
} |
810 |
gpertea |
137 |
while (isalpha(*p)) { |
811 |
gpertea |
70 |
p++; |
812 |
gpertea |
137 |
mismatches.push_back(bi); |
813 |
gpertea |
70 |
bi++; |
814 |
|
|
} |
815 |
|
|
if (*p=='^') { //reference deletion |
816 |
|
|
p++; |
817 |
gpertea |
137 |
while (isalpha(*p)) { //insert read bases |
818 |
gpertea |
70 |
p++; bi++; |
819 |
gpertea |
137 |
} |
820 |
gpertea |
70 |
} |
821 |
|
|
} |
822 |
|
|
} |
823 |
|
|
else |
824 |
|
|
{ |
825 |
|
|
//fprintf(stderr, "%s attribute not supported\n", tuple_fields[0].c_str()); |
826 |
|
|
//return false; |
827 |
|
|
} |
828 |
|
|
} |
829 |
|
|
} |
830 |
|
|
/* By convention,the NM field of the SAM record |
831 |
|
|
* counts an insertion or deletion. I dont' think |
832 |
|
|
* we want the mismatch count in the BowtieHit |
833 |
|
|
* record to reflect this. Therefore, subtract out |
834 |
|
|
* the mismatches due to in/dels |
835 |
|
|
*/ |
836 |
|
|
for(vector<CigarOp>::const_iterator itr = cigar.begin(); itr != cigar.end(); ++itr){ |
837 |
gpertea |
137 |
switch (itr->opcode) |
838 |
|
|
{ |
839 |
|
|
case MATCH: |
840 |
|
|
case REF_SKIP: |
841 |
|
|
case PAD: |
842 |
|
|
gspan += itr->length; |
843 |
|
|
break; |
844 |
|
|
case DEL: |
845 |
|
|
gspan += itr->length; |
846 |
|
|
sam_nm -= itr->length; |
847 |
|
|
break; |
848 |
|
|
case INS: |
849 |
|
|
sam_nm -= itr->length; |
850 |
|
|
break; |
851 |
gpertea |
70 |
} |
852 |
gpertea |
137 |
} |
853 |
|
|
return gspan; |
854 |
gpertea |
70 |
} |
855 |
|
|
|
856 |
gpertea |
29 |
bool SAMHitFactory::get_hit_from_buf(const char* orig_bwt_buf, |
857 |
|
|
BowtieHit& bh, |
858 |
|
|
bool strip_slash, |
859 |
|
|
char* name_out, |
860 |
|
|
char* name_tags, |
861 |
|
|
char* seq, |
862 |
|
|
char* qual) |
863 |
|
|
{ |
864 |
|
|
if (!orig_bwt_buf || !*orig_bwt_buf) |
865 |
|
|
return false; |
866 |
|
|
char bwt_buf[2048]; |
867 |
|
|
strcpy(bwt_buf, orig_bwt_buf); |
868 |
|
|
// Are we still in the header region? |
869 |
|
|
if (bwt_buf[0] == '@') |
870 |
|
|
return false; |
871 |
|
|
|
872 |
|
|
char* buf = bwt_buf; |
873 |
|
|
char* name = get_token((char**)&buf,"\t"); |
874 |
|
|
char* sam_flag_str = get_token((char**)&buf,"\t"); |
875 |
|
|
char* text_name = get_token((char**)&buf,"\t"); |
876 |
|
|
char* text_offset_str = get_token((char**)&buf,"\t"); |
877 |
|
|
const char* map_qual_str = get_token((char**)&buf,"\t"); |
878 |
|
|
char* cigar_str = get_token((char**)&buf,"\t"); |
879 |
|
|
const char* mate_ref_str = get_token((char**)&buf,"\t"); |
880 |
|
|
const char* mate_pos_str = get_token((char**)&buf,"\t"); |
881 |
|
|
const char* inferred_insert_sz_str = get_token((char**)&buf,"\t"); |
882 |
|
|
|
883 |
|
|
const char* seq_str = get_token((char**)&buf,"\t"); |
884 |
|
|
if (seq) |
885 |
|
|
strcpy(seq, seq_str); |
886 |
|
|
|
887 |
|
|
const char* qual_str = get_token((char**)&buf,"\t"); |
888 |
|
|
if (qual) |
889 |
|
|
strcpy(qual, qual_str); |
890 |
|
|
|
891 |
|
|
if (!name || |
892 |
|
|
!sam_flag_str || |
893 |
|
|
!text_name || |
894 |
|
|
!text_offset_str || |
895 |
|
|
!map_qual_str || |
896 |
|
|
!cigar_str || |
897 |
|
|
!mate_ref_str || |
898 |
|
|
!mate_pos_str || |
899 |
|
|
!inferred_insert_sz_str || |
900 |
|
|
!seq_str || |
901 |
|
|
!qual_str) |
902 |
|
|
{ |
903 |
|
|
// truncated or malformed SAM record |
904 |
|
|
return false; |
905 |
|
|
} |
906 |
|
|
|
907 |
|
|
|
908 |
|
|
int sam_flag = atoi(sam_flag_str); |
909 |
|
|
int text_offset = atoi(text_offset_str); |
910 |
|
|
|
911 |
|
|
bool end = true; |
912 |
|
|
unsigned int seg_offset = 0; |
913 |
|
|
unsigned int seg_num = 0; |
914 |
|
|
unsigned int num_segs = 0; |
915 |
|
|
|
916 |
|
|
// Copy the tag out of the name field before we might wipe it out |
917 |
gpertea |
70 |
parseSegReadName(name, name_tags, strip_slash, end, seg_offset, seg_num, num_segs); |
918 |
gpertea |
29 |
|
919 |
|
|
vector<CigarOp> cigar; |
920 |
|
|
bool spliced_alignment = false; |
921 |
gpertea |
137 |
int refspan=parseCigar(cigar, cigar_str, spliced_alignment); |
922 |
|
|
if (refspan==0) |
923 |
gpertea |
70 |
return false; |
924 |
gpertea |
29 |
//vector<string> attributes; |
925 |
|
|
//tokenize(tag_buf, " \t",attributes); |
926 |
gpertea |
70 |
|
927 |
gpertea |
29 |
bool antisense_splice = false; |
928 |
gpertea |
70 |
int sam_nm = 0; //the value of the NM tag (edit distance) |
929 |
gpertea |
137 |
//int mismatches[1024];//array with mismatch positions on the read (0-based from the left aligned end of the read) |
930 |
|
|
vector<int> mismatches; |
931 |
|
|
getSAMmismatches(buf, cigar, mismatches, sam_nm, antisense_splice); |
932 |
gpertea |
29 |
|
933 |
|
|
if (spliced_alignment) |
934 |
|
|
{ |
935 |
|
|
bh = create_hit(name, |
936 |
|
|
text_name, |
937 |
|
|
text_offset - 1, |
938 |
|
|
cigar, |
939 |
|
|
sam_flag & 0x0010, |
940 |
|
|
antisense_splice, |
941 |
gpertea |
137 |
mismatches.size(), |
942 |
gpertea |
70 |
0, |
943 |
gpertea |
29 |
end); |
944 |
|
|
} |
945 |
|
|
else |
946 |
|
|
{ |
947 |
|
|
//assert(cigar.size() == 1 && cigar[0].opcode == MATCH); |
948 |
|
|
bh = create_hit(name, |
949 |
|
|
text_name, |
950 |
|
|
text_offset - 1, // SAM files are 1-indexed |
951 |
|
|
cigar, |
952 |
|
|
sam_flag & 0x0010, |
953 |
|
|
false, |
954 |
gpertea |
137 |
mismatches.size(), |
955 |
gpertea |
29 |
0, |
956 |
|
|
end); |
957 |
|
|
} |
958 |
gpertea |
70 |
return true; |
959 |
gpertea |
29 |
} |
960 |
|
|
|
961 |
gpertea |
137 |
void cigar_add(vector<CigarOp>& cigar, CigarOp& op) { |
962 |
|
|
if (op.length<=0) return; |
963 |
|
|
if (cigar.size()>0 && cigar.back().opcode==op.opcode) { |
964 |
|
|
cigar.back().length+=op.length; |
965 |
|
|
} |
966 |
|
|
cigar.push_back(op); |
967 |
|
|
} |
968 |
gpertea |
70 |
|
969 |
gpertea |
137 |
int spliceCigar(vector<CigarOp>& splcigar, vector<CigarOp>& cigar, vector<int> mismatches, |
970 |
|
|
int &left, int spl_start, int spl_len, CigarOpCode spl_code) { |
971 |
gpertea |
70 |
//merge the original 'cigar' with the new insert/gap operation |
972 |
gpertea |
137 |
//at position spl_start and place the result into splcigar; |
973 |
|
|
int num_mismatches=mismatches.size(); |
974 |
|
|
//these offsets are relative to the beginning of alignment on reference |
975 |
|
|
int spl_ofs=spl_start-left; //relative position of splice op |
976 |
|
|
int spl_ofs_end=spl_ofs; //relative position of first ref base AFTER splice op |
977 |
|
|
CigarOp gapop(spl_code, spl_len); //for DEL or REF_SKIP |
978 |
|
|
if (spl_code==INS) |
979 |
|
|
spl_ofs_end += spl_len; |
980 |
|
|
int ref_ofs=0; //working offset on reference |
981 |
|
|
int read_ofs=0; //working offset on the read, relative to the leftmost aligned base |
982 |
|
|
bool xfound=false; |
983 |
|
|
//if (left<=spl_start+spl_len) { |
984 |
|
|
if (spl_ofs_end>0) { |
985 |
|
|
int prev_opcode=0; |
986 |
|
|
int prev_oplen=0; |
987 |
|
|
for (size_t c = 0 ; c < cigar.size(); ++c) |
988 |
|
|
{ |
989 |
|
|
int prev_read_ofs=read_ofs; |
990 |
|
|
int cur_op_ofs=ref_ofs; |
991 |
|
|
int cur_opcode=cigar[c].opcode; |
992 |
|
|
int cur_oplen=cigar[c].length; |
993 |
|
|
switch (cur_opcode) { |
994 |
|
|
case MATCH: |
995 |
|
|
ref_ofs+=cur_oplen; |
996 |
|
|
read_ofs+=cur_oplen; |
997 |
|
|
break; |
998 |
|
|
case DEL: |
999 |
|
|
case REF_SKIP: |
1000 |
|
|
case PAD: |
1001 |
|
|
ref_ofs+=cur_oplen; |
1002 |
|
|
break; |
1003 |
|
|
case SOFT_CLIP: |
1004 |
|
|
case INS: |
1005 |
|
|
read_ofs+=cur_oplen; |
1006 |
|
|
break; |
1007 |
|
|
//case HARD_CLIP: |
1008 |
|
|
} |
1009 |
|
|
if (cur_op_ofs>=spl_ofs_end || ref_ofs<=spl_ofs) { |
1010 |
|
|
if (spl_code!=INS && cur_op_ofs==spl_ofs_end) |
1011 |
|
|
{ |
1012 |
|
|
//we have to insert the gap here first |
1013 |
|
|
cigar_add(splcigar, gapop); |
1014 |
|
|
} |
1015 |
|
|
cigar_add(splcigar, cigar[c]); |
1016 |
|
|
} |
1017 |
|
|
else //if (ref_ofs>spl_ofs) { |
1018 |
|
|
{ //op intersection |
1019 |
|
|
xfound=true; |
1020 |
|
|
if (spl_code==INS) { |
1021 |
|
|
//we have to shorten cur_opcode |
1022 |
|
|
// find the overlap between current range |
1023 |
|
|
int ovl_start = (cur_op_ofs>spl_ofs) ? cur_op_ofs : spl_ofs; |
1024 |
|
|
int ovl_end = (ref_ofs>spl_ofs_end) ? spl_ofs_end : ref_ofs; |
1025 |
|
|
CigarOp op(cigar[c]); |
1026 |
|
|
int len=op.length; |
1027 |
|
|
len-=(ovl_end-ovl_start); |
1028 |
|
|
if (len>0) { |
1029 |
|
|
op.length=len; |
1030 |
|
|
cigar_add(splcigar, op); |
1031 |
|
|
} |
1032 |
|
|
} |
1033 |
|
|
else {//DEL or REF_SKIP |
1034 |
|
|
//spl_ofs == spl_ofs_end |
1035 |
|
|
//we have to split cur_opcode |
1036 |
|
|
CigarOp op(cigar[c]); |
1037 |
|
|
op.length=spl_ofs-cur_op_ofs; |
1038 |
|
|
cigar_add(splcigar, op); |
1039 |
|
|
cigar_add(splcigar, gapop); |
1040 |
|
|
op.length=ref_ofs-spl_ofs; |
1041 |
|
|
cigar_add(splcigar,op); |
1042 |
|
|
} |
1043 |
|
|
} //op intersection |
1044 |
|
|
prev_opcode=cur_opcode; |
1045 |
|
|
prev_oplen=cur_oplen; |
1046 |
|
|
} //for each cigar opcode |
1047 |
|
|
} //intersection possible |
1048 |
|
|
if (!xfound) {//no intersection found between splice event and alignment |
1049 |
|
|
if (spl_ofs_end<=0) { |
1050 |
|
|
//alignment starts after the splice event |
1051 |
|
|
if (spl_code==INS) left-=spl_len; |
1052 |
|
|
else left+=spl_len; |
1053 |
|
|
} |
1054 |
|
|
//else { |
1055 |
|
|
//alignment ends before the splice event |
1056 |
|
|
//nothing to do |
1057 |
|
|
// } |
1058 |
|
|
} |
1059 |
gpertea |
97 |
|
1060 |
gpertea |
137 |
//TODO: adjust mismatches if opcode==INS and they fall in the removed part |
1061 |
|
|
//TODO: what about anchor mismatches for REF_SKIP |
1062 |
|
|
return num_mismatches; |
1063 |
gpertea |
70 |
} |
1064 |
|
|
|
1065 |
gpertea |
69 |
bool SplicedSAMHitFactory::get_hit_from_buf(const char* orig_bwt_buf, |
1066 |
|
|
BowtieHit& bh, |
1067 |
|
|
bool strip_slash, |
1068 |
|
|
char* name_out, |
1069 |
|
|
char* name_tags, |
1070 |
|
|
char* seq, |
1071 |
|
|
char* qual) |
1072 |
|
|
{ |
1073 |
|
|
if (!orig_bwt_buf || !*orig_bwt_buf) |
1074 |
|
|
return false; |
1075 |
gpertea |
29 |
|
1076 |
gpertea |
69 |
char bwt_buf[2048]; |
1077 |
|
|
strcpy(bwt_buf, orig_bwt_buf); |
1078 |
|
|
// Are we still in the header region? |
1079 |
|
|
if (bwt_buf[0] == '@') |
1080 |
|
|
return false; |
1081 |
gpertea |
29 |
|
1082 |
gpertea |
69 |
char* buf = bwt_buf; |
1083 |
|
|
char* name = get_token((char**)&buf,"\t"); |
1084 |
|
|
char* sam_flag_str = get_token((char**)&buf,"\t"); |
1085 |
|
|
char* text_name = get_token((char**)&buf,"\t"); |
1086 |
|
|
char* text_offset_str = get_token((char**)&buf,"\t"); |
1087 |
|
|
const char* map_qual_str = get_token((char**)&buf,"\t"); |
1088 |
|
|
char* cigar_str = get_token((char**)&buf,"\t"); |
1089 |
|
|
const char* mate_ref_str = get_token((char**)&buf,"\t"); |
1090 |
|
|
const char* mate_pos_str = get_token((char**)&buf,"\t"); |
1091 |
|
|
const char* inferred_insert_sz_str = get_token((char**)&buf,"\t"); |
1092 |
gpertea |
137 |
//int num_mismatches=0; |
1093 |
|
|
//int mismatches[1024]; //list of 0-based mismatch positions in this read |
1094 |
gpertea |
70 |
//parsed from SAM's MD:Z: tag |
1095 |
gpertea |
69 |
const char* seq_str = get_token((char**)&buf,"\t"); |
1096 |
|
|
if (seq) |
1097 |
|
|
strcpy(seq, seq_str); |
1098 |
|
|
|
1099 |
|
|
const char* qual_str = get_token((char**)&buf,"\t"); |
1100 |
|
|
if (qual) |
1101 |
|
|
strcpy(qual, qual_str); |
1102 |
|
|
|
1103 |
|
|
if (!name || |
1104 |
|
|
!sam_flag_str || |
1105 |
|
|
!text_name || |
1106 |
|
|
!text_offset_str || |
1107 |
|
|
!map_qual_str || |
1108 |
|
|
!cigar_str || |
1109 |
|
|
!mate_ref_str || |
1110 |
|
|
!mate_pos_str || |
1111 |
|
|
!inferred_insert_sz_str || |
1112 |
|
|
!seq_str || |
1113 |
|
|
!qual_str) |
1114 |
|
|
{ |
1115 |
|
|
// truncated or malformed SAM record |
1116 |
|
|
return false; |
1117 |
|
|
} |
1118 |
|
|
|
1119 |
|
|
|
1120 |
|
|
int sam_flag = atoi(sam_flag_str); |
1121 |
|
|
int text_offset = atoi(text_offset_str); |
1122 |
gpertea |
137 |
text_offset--; //make it 0-based (SAM is 1-based, Bowtie is 0-based) |
1123 |
gpertea |
69 |
bool end = true; |
1124 |
|
|
unsigned int seg_offset = 0; |
1125 |
|
|
unsigned int seg_num = 0; |
1126 |
|
|
unsigned int num_segs = 0; |
1127 |
|
|
|
1128 |
|
|
// Copy the tag out of the name field before we might wipe it out |
1129 |
gpertea |
70 |
parseSegReadName(name, name_tags, strip_slash, end, seg_offset, seg_num, num_segs); |
1130 |
gpertea |
69 |
|
1131 |
gpertea |
70 |
vector<CigarOp> samcigar; |
1132 |
|
|
bool spliced_alignment = false; |
1133 |
gpertea |
137 |
int refspan=parseCigar(samcigar, cigar_str, spliced_alignment); |
1134 |
|
|
if (refspan==0) |
1135 |
|
|
return false; |
1136 |
gpertea |
70 |
bool antisense_splice = false; |
1137 |
|
|
int sam_nm = 0; |
1138 |
gpertea |
137 |
vector<int> mismatches; |
1139 |
|
|
getSAMmismatches(buf, samcigar, mismatches, sam_nm, antisense_splice); |
1140 |
gpertea |
69 |
|
1141 |
gpertea |
70 |
//############################################## |
1142 |
gpertea |
69 |
|
1143 |
|
|
// Add this alignment to the table of hits for this half of the |
1144 |
|
|
// Bowtie map |
1145 |
|
|
|
1146 |
|
|
// Parse the text_name field to recover the splice coords |
1147 |
|
|
vector<string> toks; |
1148 |
|
|
tokenize_strict(text_name, "|", toks); |
1149 |
|
|
|
1150 |
|
|
int num_extra_toks = (int)toks.size() - 6; |
1151 |
|
|
|
1152 |
|
|
if (num_extra_toks >= 0) |
1153 |
|
|
{ |
1154 |
|
|
static const uint8_t left_window_edge_field = 1; |
1155 |
|
|
static const uint8_t splice_field = 2; |
1156 |
|
|
//static const uint8_t right_window_edge_field = 3; |
1157 |
|
|
static const uint8_t junction_type_field = 4; |
1158 |
|
|
static const uint8_t strand_field = 5; |
1159 |
|
|
|
1160 |
|
|
string contig = toks[0]; |
1161 |
|
|
for (int t = 1; t <= num_extra_toks; ++t) |
1162 |
|
|
{ |
1163 |
|
|
contig += "|"; |
1164 |
|
|
contig += toks[t]; |
1165 |
|
|
} |
1166 |
|
|
|
1167 |
|
|
vector<string> splice_toks; |
1168 |
|
|
tokenize(toks[num_extra_toks + splice_field], "-", splice_toks); |
1169 |
|
|
if (splice_toks.size() != 2) |
1170 |
|
|
{ |
1171 |
|
|
fprintf(stderr, "Warning: found malformed splice record, skipping:\n"); |
1172 |
gpertea |
70 |
//fprintf(stderr, "\t%s (token: %s)\n", text_name, |
1173 |
gpertea |
69 |
// toks[num_extra_toks + splice_field].c_str()); |
1174 |
|
|
return false; |
1175 |
|
|
} |
1176 |
gpertea |
137 |
|
1177 |
|
|
string junction_strand = toks[num_extra_toks + strand_field]; |
1178 |
|
|
if(junction_strand != "rev" && junction_strand != "fwd"){ |
1179 |
|
|
fprintf(stderr, "Malformed insertion record\n"); |
1180 |
|
|
return false; |
1181 |
|
|
} |
1182 |
|
|
|
1183 |
|
|
int spl_num_mismatches=0; |
1184 |
|
|
|
1185 |
gpertea |
69 |
// |
1186 |
|
|
// check for an insertion hit |
1187 |
|
|
// |
1188 |
|
|
if(toks[num_extra_toks + junction_type_field] == "ins") |
1189 |
|
|
{ |
1190 |
gpertea |
137 |
//int8_t spliced_read_len = strlen(seq_str); |
1191 |
|
|
//TODO FIXME: use the CIGAR instead of seq length! |
1192 |
gpertea |
70 |
// The 0-based position of the left edge of the alignment. Note that this |
1193 |
gpertea |
137 |
// value may need to be further corrected to account for the presence of |
1194 |
gpertea |
70 |
// of the insertion. |
1195 |
gpertea |
137 |
int left = atoi(toks[num_extra_toks + left_window_edge_field].c_str()) + text_offset; |
1196 |
gpertea |
69 |
|
1197 |
gpertea |
70 |
// The 0-based position of the last genomic sequence before the insertion |
1198 |
gpertea |
137 |
int left_splice_pos = atoi(splice_toks[0].c_str()); |
1199 |
gpertea |
69 |
|
1200 |
|
|
string insertedSequence = splice_toks[1]; |
1201 |
gpertea |
70 |
// The 0-based position of the first genomic sequence after the insertion |
1202 |
|
|
|
1203 |
gpertea |
137 |
vector<CigarOp> splcigar; |
1204 |
|
|
|
1205 |
|
|
//TODO: implement this |
1206 |
|
|
//this also updates left to the adjusted genomic coordinates |
1207 |
|
|
spl_num_mismatches=spliceCigar(splcigar, samcigar, mismatches, |
1208 |
|
|
left, left_splice_pos, insertedSequence.length(), INS); |
1209 |
|
|
/* |
1210 |
gpertea |
69 |
uint32_t right_splice_pos = left_splice_pos + 1; |
1211 |
gpertea |
137 |
|
1212 |
|
|
//uint32_t right = left + spliced_read_len - 1; |
1213 |
|
|
int right = left + refspan - 1; |
1214 |
|
|
|
1215 |
gpertea |
69 |
if(left > left_splice_pos){ |
1216 |
gpertea |
137 |
//The genomic position of the left edge of the alignment needs to be corrected |
1217 |
|
|
//If the alignment does not extend into the insertion, simply subtract the length |
1218 |
|
|
//of the inserted sequence, otherwise, just set it equal to the right edge |
1219 |
gpertea |
69 |
left = left - insertedSequence.length(); |
1220 |
|
|
if(left < right_splice_pos){ |
1221 |
|
|
left = right_splice_pos; |
1222 |
|
|
} |
1223 |
|
|
} |
1224 |
|
|
if(right > left_splice_pos){ |
1225 |
|
|
right = right - insertedSequence.length(); |
1226 |
|
|
if(right < left_splice_pos){ |
1227 |
|
|
right = left_splice_pos; |
1228 |
|
|
} |
1229 |
|
|
} |
1230 |
gpertea |
137 |
// Now, right and left should be properly transformed into genomic coordinates |
1231 |
|
|
// We should be able to deduce how much the alignment matches the insertion |
1232 |
|
|
// simply based on the length of the read |
1233 |
gpertea |
69 |
int left_match_length = 0; |
1234 |
|
|
if(left <= left_splice_pos){ |
1235 |
|
|
left_match_length = left_splice_pos - left + 1; |
1236 |
|
|
} |
1237 |
|
|
int right_match_length = 0; |
1238 |
|
|
if(right >= right_splice_pos){ |
1239 |
|
|
right_match_length = right - right_splice_pos + 1; |
1240 |
|
|
} |
1241 |
|
|
int insertion_match_length = spliced_read_len - left_match_length - right_match_length; |
1242 |
|
|
|
1243 |
|
|
if(left_match_length <= 0 || right_match_length <= 0 || insertion_match_length <= 0) |
1244 |
|
|
return false; |
1245 |
|
|
|
1246 |
gpertea |
70 |
//char* pch = strtok( mismatches, ","); |
1247 |
|
|
//unsigned char num_mismatches = 0; |
1248 |
|
|
//text_offset holds the left end of the alignment, |
1249 |
|
|
//RELATIVE TO the start of the contig |
1250 |
gpertea |
69 |
|
1251 |
gpertea |
70 |
//The 0-based relative position of the left-most character |
1252 |
|
|
//before the insertion in the contig |
1253 |
gpertea |
69 |
int relative_splice_pos = left_splice_pos - atoi(toks[num_extra_toks + left_window_edge_field].c_str()); |
1254 |
gpertea |
137 |
for (size_t i=0;i<mismatches.size();++i) { |
1255 |
gpertea |
70 |
int mismatch_pos = mismatches[i]; |
1256 |
|
|
// for reversely mapped reads, |
1257 |
|
|
//find the correct mismatched position. |
1258 |
|
|
if (sam_flag & 0x0010){ |
1259 |
|
|
mismatch_pos = spliced_read_len - mismatch_pos - 1; |
1260 |
|
|
} |
1261 |
gpertea |
69 |
|
1262 |
gpertea |
70 |
//Only count mismatches outside of the insertion region |
1263 |
|
|
// If there is a mismatch within the insertion, |
1264 |
|
|
// disallow this hit |
1265 |
gpertea |
137 |
if(mismatch_pos + text_offset <= relative_splice_pos || |
1266 |
|
|
mismatch_pos + text_offset > relative_splice_pos + insertedSequence.length()){ |
1267 |
gpertea |
70 |
spl_num_mismatches++; |
1268 |
|
|
}else{ |
1269 |
|
|
return false; |
1270 |
gpertea |
69 |
} |
1271 |
|
|
} |
1272 |
gpertea |
137 |
*/ |
1273 |
|
|
//vector<CigarOp> splcigar; |
1274 |
|
|
//spliceCigar(splcigar, samcigar, left_match_length, insertion_match_length, right_match_length, INS); |
1275 |
|
|
//splcigar.push_back(CigarOp(MATCH, left_match_length)); |
1276 |
|
|
//splcigar.push_back(CigarOp(INS, insertion_match_length)); |
1277 |
|
|
//splcigar.push_back(CigarOp(MATCH, right_match_length)); |
1278 |
gpertea |
69 |
|
1279 |
|
|
bh = create_hit(name, |
1280 |
|
|
contig, |
1281 |
|
|
left, |
1282 |
gpertea |
70 |
//splcigar, |
1283 |
|
|
splcigar, |
1284 |
|
|
sam_flag & 0x0010, |
1285 |
gpertea |
69 |
junction_strand == "rev", |
1286 |
gpertea |
70 |
spl_num_mismatches, |
1287 |
gpertea |
69 |
0, |
1288 |
|
|
end); |
1289 |
|
|
return true; |
1290 |
gpertea |
70 |
} //"ins" |
1291 |
gpertea |
137 |
else //"del" or intron |
1292 |
gpertea |
69 |
{ |
1293 |
gpertea |
137 |
int left = atoi(toks[num_extra_toks + left_window_edge_field].c_str()) + text_offset; |
1294 |
|
|
int spliced_read_len = strlen(seq_str); //FIXME TODO: use CIGAR instead |
1295 |
|
|
int left_splice_pos = atoi(splice_toks[0].c_str()) - left + 1; |
1296 |
gpertea |
69 |
if(left_splice_pos > spliced_read_len) left_splice_pos = spliced_read_len; |
1297 |
gpertea |
137 |
int right_splice_pos = spliced_read_len - left_splice_pos; |
1298 |
gpertea |
69 |
|
1299 |
|
|
if (right_splice_pos <= 0 || left_splice_pos <= 0) |
1300 |
|
|
return false; |
1301 |
|
|
|
1302 |
|
|
if ((sam_flag & 0x0010) == 0) //###### |
1303 |
|
|
{ |
1304 |
gpertea |
135 |
if (left_splice_pos + seg_offset < _anchor_length) |
1305 |
|
|
return false; |
1306 |
|
|
} |
1307 |
|
|
else |
1308 |
|
|
{ |
1309 |
|
|
if (right_splice_pos + seg_offset < _anchor_length) |
1310 |
|
|
return false; |
1311 |
|
|
} |
1312 |
gpertea |
69 |
//uint32_t right = atoi(splice_toks[1].c_str()) + right_splice_pos; |
1313 |
|
|
//atoi(toks[right_window_edge_field].c_str()); |
1314 |
|
|
int gap_len = atoi(splice_toks[1].c_str()) - atoi(splice_toks[0].c_str()) - 1; |
1315 |
|
|
|
1316 |
|
|
int mismatches_in_anchor = 0; |
1317 |
gpertea |
137 |
for (size_t i=0;i<mismatches.size();++i) { |
1318 |
|
|
spl_num_mismatches++; |
1319 |
|
|
int mismatch_pos = mismatches[i]; |
1320 |
|
|
if (((sam_flag & 0x0010) == 0 && abs(mismatch_pos - left_splice_pos) < (int)min_anchor_len) || |
1321 |
|
|
((sam_flag & 0x0010) != 0 && |
1322 |
gpertea |
70 |
abs(((int)spliced_read_len - left_splice_pos + 1) - mismatch_pos)) < (int)min_anchor_len) |
1323 |
|
|
mismatches_in_anchor++; |
1324 |
gpertea |
137 |
} |
1325 |
gpertea |
69 |
|
1326 |
|
|
// FIXME: we probably should exclude these hits somewhere, but this |
1327 |
|
|
// isn't the right place |
1328 |
gpertea |
70 |
vector<CigarOp> splcigar; |
1329 |
|
|
CigarOpCode opcode=(toks[num_extra_toks + junction_type_field] == "del")? DEL : REF_SKIP ; |
1330 |
gpertea |
137 |
//TODO: FIXME: use relative offsets for this call instead, like for the INS case |
1331 |
|
|
spl_num_mismatches=spliceCigar(splcigar, samcigar, mismatches, left, left_splice_pos, gap_len, opcode); |
1332 |
gpertea |
70 |
/* |
1333 |
|
|
splcigar.push_back(CigarOp(MATCH, left_splice_pos)); |
1334 |
gpertea |
69 |
if(toks[num_extra_toks + junction_type_field] == "del"){ |
1335 |
gpertea |
70 |
splcigar.push_back(CigarOp(DEL, gap_len)); |
1336 |
gpertea |
69 |
}else{ |
1337 |
gpertea |
70 |
splcigar.push_back(CigarOp(REF_SKIP, gap_len)); |
1338 |
gpertea |
69 |
} |
1339 |
gpertea |
70 |
splcigar.push_back(CigarOp(MATCH, right_splice_pos)); |
1340 |
|
|
*/ |
1341 |
gpertea |
69 |
bh = create_hit(name, |
1342 |
|
|
contig, |
1343 |
|
|
left, |
1344 |
gpertea |
70 |
splcigar, |
1345 |
|
|
(sam_flag & 0x0010), |
1346 |
gpertea |
69 |
junction_strand == "rev", |
1347 |
gpertea |
70 |
spl_num_mismatches, |
1348 |
gpertea |
69 |
mismatches_in_anchor, |
1349 |
|
|
end); |
1350 |
|
|
return true; |
1351 |
|
|
} |
1352 |
gpertea |
70 |
} //parse splice data |
1353 |
gpertea |
69 |
else |
1354 |
|
|
{ |
1355 |
|
|
fprintf(stderr, "Warning: found malformed splice record, skipping\n"); |
1356 |
|
|
//fprintf(stderr, "%s\n", orig_bwt_buf); |
1357 |
|
|
// continue; |
1358 |
|
|
return false; |
1359 |
|
|
} |
1360 |
|
|
|
1361 |
|
|
return false; |
1362 |
|
|
} |
1363 |
|
|
|
1364 |
|
|
|
1365 |
gpertea |
135 |
|
1366 |
gpertea |
29 |
bool BAMHitFactory::get_hit_from_buf(const char* orig_bwt_buf, |
1367 |
|
|
BowtieHit& bh, bool strip_slash, |
1368 |
|
|
char* name_out, char* name_tags, |
1369 |
|
|
char* seq, char* qual) { |
1370 |
|
|
if (_sam_header==NULL) |
1371 |
|
|
err_die("Error: no SAM header when BAMHitFactory::get_hit_from_buf()!"); |
1372 |
|
|
const bam1_t* hit_buf = (const bam1_t*)orig_bwt_buf; |
1373 |
|
|
|
1374 |
|
|
|
1375 |
|
|
uint32_t sam_flag = hit_buf->core.flag; |
1376 |
|
|
|
1377 |
|
|
int text_offset = hit_buf->core.pos; |
1378 |
|
|
int text_mate_pos = hit_buf->core.mpos; |
1379 |
|
|
int target_id = hit_buf->core.tid; |
1380 |
|
|
int mate_target_id = hit_buf->core.mtid; |
1381 |
|
|
|
1382 |
|
|
vector<CigarOp> cigar; |
1383 |
|
|
bool spliced_alignment = false; |
1384 |
|
|
int num_hits = 1; |
1385 |
|
|
|
1386 |
|
|
double mapQ = hit_buf->core.qual; |
1387 |
|
|
long double error_prob; |
1388 |
|
|
if (mapQ > 0) |
1389 |
|
|
{ |
1390 |
|
|
long double p = (-1.0 * mapQ) / 10.0; |
1391 |
|
|
error_prob = pow(10.0L, p); |
1392 |
|
|
} |
1393 |
|
|
else |
1394 |
|
|
{ |
1395 |
|
|
error_prob = 1.0; |
1396 |
|
|
} |
1397 |
|
|
|
1398 |
|
|
//header->target_name[c->tid] |
1399 |
|
|
|
1400 |
|
|
bool end = true; |
1401 |
|
|
unsigned int seg_offset = 0; |
1402 |
|
|
unsigned int seg_num = 0; |
1403 |
|
|
unsigned int num_segs = 0; |
1404 |
|
|
// Copy the tag out of the name field before we might wipe it out |
1405 |
|
|
char* qname = bam1_qname(hit_buf); |
1406 |
|
|
char* pipe = strrchr(qname, '|'); |
1407 |
|
|
if (pipe) |
1408 |
|
|
{ |
1409 |
|
|
if (name_tags) |
1410 |
|
|
strcpy(name_tags, pipe); |
1411 |
|
|
|
1412 |
|
|
char* tag_buf = pipe + 1; |
1413 |
|
|
if (strchr(tag_buf, ':')) |
1414 |
|
|
{ |
1415 |
|
|
sscanf(tag_buf, "%u:%u:%u", &seg_offset, &seg_num, &num_segs); |
1416 |
|
|
if (seg_num + 1 == num_segs) |
1417 |
|
|
end = true; |
1418 |
|
|
else |
1419 |
|
|
end = false; |
1420 |
|
|
} |
1421 |
|
|
|
1422 |
|
|
*pipe = 0; |
1423 |
|
|
} |
1424 |
|
|
|
1425 |
|
|
|
1426 |
|
|
if (target_id < 0) { |
1427 |
|
|
//assert(cigar.size() == 1 && cigar[0].opcode == MATCH); |
1428 |
|
|
bh = create_hit(qname, |
1429 |
|
|
"*", //ref_name |
1430 |
|
|
0, //left coord |
1431 |
|
|
0, //read_len |
1432 |
|
|
false, //antisense_aln |
1433 |
|
|
0, //edit_dist |
1434 |
|
|
end); |
1435 |
|
|
return true; |
1436 |
|
|
} |
1437 |
|
|
|
1438 |
|
|
//string text_name = ((samfile_t*)(hs._hit_file))->header->target_name[target_id]; |
1439 |
|
|
string text_name = _sam_header->target_name[target_id]; |
1440 |
|
|
for (int i = 0; i < hit_buf->core.n_cigar; ++i) |
1441 |
|
|
{ |
1442 |
|
|
//char* t; |
1443 |
|
|
|
1444 |
|
|
int length = bam1_cigar(hit_buf)[i] >> BAM_CIGAR_SHIFT; |
1445 |
|
|
if (length <= 0) |
1446 |
|
|
{ |
1447 |
|
|
fprintf (stderr, "BAM error: CIGAR op has zero length\n"); |
1448 |
|
|
return false; |
1449 |
|
|
} |
1450 |
|
|
|
1451 |
|
|
CigarOpCode opcode; |
1452 |
|
|
switch(bam1_cigar(hit_buf)[i] & BAM_CIGAR_MASK) |
1453 |
|
|
{ |
1454 |
|
|
case BAM_CMATCH: opcode = MATCH; break; |
1455 |
|
|
case BAM_CINS: opcode = INS; break; |
1456 |
|
|
case BAM_CDEL: opcode = DEL; break; |
1457 |
|
|
case BAM_CSOFT_CLIP: opcode = SOFT_CLIP; break; |
1458 |
|
|
case BAM_CHARD_CLIP: opcode = HARD_CLIP; break; |
1459 |
|
|
case BAM_CPAD: opcode = PAD; break; |
1460 |
|
|
case BAM_CREF_SKIP: |
1461 |
|
|
opcode = REF_SKIP; |
1462 |
|
|
spliced_alignment = true; |
1463 |
|
|
if (length > (int)max_report_intron_length) |
1464 |
|
|
{ |
1465 |
|
|
//fprintf(stderr, "Encounter REF_SKIP > max_gene_length, skipping\n"); |
1466 |
|
|
return false; |
1467 |
|
|
} |
1468 |
|
|
break; |
1469 |
|
|
default: |
1470 |
|
|
fprintf (stderr, "BAM read: invalid CIGAR operation\n"); |
1471 |
|
|
return false; |
1472 |
|
|
} |
1473 |
|
|
if (opcode != HARD_CLIP) |
1474 |
|
|
cigar.push_back(CigarOp(opcode, length)); |
1475 |
|
|
} |
1476 |
|
|
|
1477 |
|
|
string mrnm; |
1478 |
|
|
if (mate_target_id >= 0) { |
1479 |
|
|
if (mate_target_id == target_id) { |
1480 |
|
|
//mrnm = ((samfile_t*)(hs._hit_file))->header->target_name[mate_target_id]; |
1481 |
|
|
mrnm = _sam_header->target_name[mate_target_id]; |
1482 |
|
|
} |
1483 |
|
|
else { |
1484 |
|
|
//fprintf(stderr, "Trans-spliced mates are not currently supported, skipping\n"); |
1485 |
|
|
return false; |
1486 |
|
|
} |
1487 |
|
|
} |
1488 |
|
|
else { |
1489 |
|
|
text_mate_pos = 0; |
1490 |
|
|
} |
1491 |
|
|
//CuffStrand source_strand = CUFF_STRAND_UNKNOWN; |
1492 |
|
|
bool antisense_splice=false; |
1493 |
|
|
unsigned char num_mismatches = 0; |
1494 |
|
|
unsigned char num_splice_anchor_mismatches = 0; |
1495 |
|
|
|
1496 |
|
|
uint8_t* ptr = bam_aux_get(hit_buf, "XS"); |
1497 |
|
|
if (ptr) { |
1498 |
|
|
char src_strand_char = bam_aux2A(ptr); |
1499 |
|
|
if (src_strand_char == '-') |
1500 |
|
|
antisense_splice = true; |
1501 |
|
|
// else if (src_strand_char == '+') |
1502 |
|
|
// source_strand = CUFF_FWD; |
1503 |
|
|
} |
1504 |
|
|
|
1505 |
|
|
ptr = bam_aux_get(hit_buf, "NM"); |
1506 |
|
|
if (ptr) { |
1507 |
|
|
num_mismatches = bam_aux2i(ptr); |
1508 |
|
|
} |
1509 |
|
|
|
1510 |
|
|
ptr = bam_aux_get(hit_buf, "NH"); |
1511 |
|
|
if (ptr) { |
1512 |
|
|
num_hits = bam_aux2i(ptr); |
1513 |
|
|
} |
1514 |
|
|
|
1515 |
|
|
//bool antisense_aln = bam1_strand(hit_buf); |
1516 |
|
|
|
1517 |
|
|
//if (_rg_props.strandedness() == STRANDED_PROTOCOL && source_strand == CUFF_STRAND_UNKNOWN) |
1518 |
|
|
// source_strand = use_stranded_protocol(sam_flag, antisense_aln, _rg_props.mate_strand_mapping()); |
1519 |
|
|
if (spliced_alignment) { |
1520 |
|
|
//if (source_strand == CUFF_STRAND_UNKNOWN) { |
1521 |
|
|
// fprintf(stderr, "BAM record error: found spliced alignment without XS attribute\n"); |
1522 |
|
|
// } |
1523 |
|
|
bh = create_hit(qname, |
1524 |
|
|
text_name, |
1525 |
|
|
text_offset, // BAM files are 0-indexed |
1526 |
|
|
cigar, |
1527 |
|
|
sam_flag & 0x0010, |
1528 |
|
|
antisense_splice, |
1529 |
|
|
num_mismatches, |
1530 |
|
|
num_splice_anchor_mismatches, |
1531 |
|
|
end); |
1532 |
|
|
|
1533 |
|
|
} |
1534 |
|
|
else { |
1535 |
|
|
//assert(_rg_props.strandedness() == STRANDED_PROTOCOL || source_strand == CUFF_STRAND_UNKNOWN); |
1536 |
|
|
//assert(cigar.size() == 1 && cigar[0].opcode == MATCH); |
1537 |
|
|
bh = create_hit(qname, |
1538 |
|
|
text_name, |
1539 |
|
|
text_offset, // BAM files are 0-indexed |
1540 |
|
|
cigar, |
1541 |
|
|
sam_flag & 0x0010, |
1542 |
|
|
false, |
1543 |
|
|
num_mismatches, |
1544 |
|
|
0, |
1545 |
|
|
end); |
1546 |
|
|
} |
1547 |
|
|
if (seq!=NULL) { |
1548 |
|
|
char *bseq = (char*)bam1_seq(hit_buf); |
1549 |
|
|
for(int i=0;i<(hit_buf->core.l_qseq);i++) { |
1550 |
|
|
char v = bam1_seqi(bseq,i); |
1551 |
|
|
seq[i]=bam_nt16_rev_table[v]; |
1552 |
|
|
} |
1553 |
|
|
seq[hit_buf->core.l_qseq]=0; |
1554 |
|
|
} |
1555 |
|
|
if (qual!=NULL) { |
1556 |
|
|
char *bq = (char*)bam1_qual(hit_buf); |
1557 |
|
|
for(int i=0;i<(hit_buf->core.l_qseq);i++) { |
1558 |
|
|
qual[i]=bq[i]+33; |
1559 |
|
|
} |
1560 |
|
|
qual[hit_buf->core.l_qseq]=0; |
1561 |
|
|
} |
1562 |
|
|
return true; |
1563 |
|
|
} |
1564 |
|
|
|
1565 |
|
|
bool BAMHitFactory::inspect_header(HitStream& hs) |
1566 |
|
|
{ |
1567 |
|
|
bam_header_t* header = ((samfile_t*)(hs._hit_file))->header; |
1568 |
|
|
|
1569 |
|
|
if (header == NULL) { |
1570 |
|
|
fprintf(stderr, "Warning: No BAM header\n"); |
1571 |
|
|
return false; |
1572 |
|
|
} |
1573 |
|
|
if (header->l_text == 0) { |
1574 |
|
|
fprintf(stderr, "Warning: BAM header has 0 length or is corrupted. Try using 'samtools reheader'.\n"); |
1575 |
|
|
return false; |
1576 |
|
|
} |
1577 |
|
|
return true; |
1578 |
|
|
} |
1579 |
|
|
|
1580 |
|
|
|
1581 |
|
|
void get_mapped_reads(FILE* bwtf, |
1582 |
|
|
HitTable& hits, |
1583 |
|
|
HitFactory& hit_factory, |
1584 |
|
|
bool strip_slash, |
1585 |
|
|
bool verbose) |
1586 |
|
|
{ |
1587 |
|
|
|
1588 |
|
|
|
1589 |
|
|
char bwt_buf[2048]; |
1590 |
|
|
uint32_t reads_extracted = 0; |
1591 |
|
|
|
1592 |
|
|
while (fgets(bwt_buf, 2048, bwtf)) |
1593 |
|
|
{ |
1594 |
|
|
// Chomp the newline |
1595 |
|
|
char* nl = strrchr(bwt_buf, '\n'); |
1596 |
|
|
if (nl) *nl = 0; |
1597 |
|
|
if (*bwt_buf == 0) |
1598 |
|
|
continue; |
1599 |
|
|
// Get a new record from the tab-delimited Bowtie map |
1600 |
|
|
BowtieHit bh; |
1601 |
|
|
if (hit_factory.get_hit_from_buf(bwt_buf, bh, strip_slash)) |
1602 |
|
|
{ |
1603 |
|
|
// Only check uniqueness if these hits are spliced |
1604 |
|
|
hits.add_hit(bh, true); |
1605 |
|
|
} |
1606 |
|
|
reads_extracted++; |
1607 |
|
|
} |
1608 |
|
|
|
1609 |
|
|
// This will sort the map by insert id. |
1610 |
|
|
hits.finalize(); |
1611 |
|
|
fprintf(stderr, "Extracted %d alignments from Bowtie map\n", reads_extracted); |
1612 |
|
|
} |
1613 |
|
|
|
1614 |
|
|
|
1615 |
|
|
/* |
1616 |
|
|
AlignStatus status(const BowtieHit* align) |
1617 |
|
|
{ |
1618 |
|
|
if (!align) |
1619 |
|
|
return UNALIGNED; |
1620 |
|
|
if (align->contiguous()) |
1621 |
|
|
return CONTIGUOUS; |
1622 |
|
|
return SPLICED; |
1623 |
|
|
} |
1624 |
|
|
*/ |
1625 |
|
|
|
1626 |
|
|
void add_hits_to_coverage(const HitList& hits, vector<unsigned short>& DoC) |
1627 |
|
|
{ |
1628 |
|
|
int max_hit_pos = -1; |
1629 |
|
|
for (size_t i = 0; i < hits.size(); ++i) |
1630 |
|
|
{ |
1631 |
|
|
max_hit_pos = max((int)hits[i].right(),max_hit_pos); |
1632 |
|
|
} |
1633 |
|
|
|
1634 |
|
|
if ((int)DoC.size() < max_hit_pos) |
1635 |
|
|
DoC.resize(max_hit_pos); |
1636 |
|
|
|
1637 |
|
|
for (size_t i = 0; i < hits.size(); ++i) |
1638 |
|
|
{ |
1639 |
|
|
const BowtieHit& bh = hits[i]; |
1640 |
|
|
|
1641 |
|
|
// split up the coverage contibution for this reads |
1642 |
|
|
size_t j = bh.left(); |
1643 |
|
|
const vector<CigarOp>& cigar = bh.cigar(); |
1644 |
|
|
|
1645 |
|
|
for (size_t c = 0 ; c < cigar.size(); ++c) |
1646 |
|
|
{ |
1647 |
|
|
switch(cigar[c].opcode) |
1648 |
|
|
{ |
1649 |
|
|
case MATCH: |
1650 |
|
|
for (size_t m = 0; m < cigar[c].length; ++m) |
1651 |
|
|
{ |
1652 |
|
|
if (DoC[j + m] < 0xFFFF) |
1653 |
|
|
DoC[j + m]++; |
1654 |
|
|
} |
1655 |
|
|
//fall through this case to REF_SKIP is intentional |
1656 |
|
|
case REF_SKIP: |
1657 |
|
|
j += cigar[c].length; |
1658 |
|
|
break; |
1659 |
|
|
default: |
1660 |
|
|
break; |
1661 |
|
|
} |
1662 |
|
|
|
1663 |
|
|
} |
1664 |
|
|
} |
1665 |
|
|
} |
1666 |
|
|
|
1667 |
|
|
void add_hit_to_coverage(const BowtieHit& bh, vector<unsigned int>& DoC) |
1668 |
|
|
{ |
1669 |
|
|
if ((int)DoC.size() < bh.right()) |
1670 |
|
|
DoC.resize(bh.right()); |
1671 |
|
|
|
1672 |
|
|
// split up the coverage contibution for this reads |
1673 |
|
|
size_t j = bh.left(); |
1674 |
|
|
const vector<CigarOp>& cigar = bh.cigar(); |
1675 |
|
|
|
1676 |
|
|
for (size_t c = 0 ; c < cigar.size(); ++c) |
1677 |
|
|
{ |
1678 |
|
|
switch(cigar[c].opcode) |
1679 |
|
|
{ |
1680 |
|
|
case MATCH: |
1681 |
|
|
for (size_t m = 0; m < cigar[c].length; ++m) |
1682 |
|
|
{ |
1683 |
gpertea |
135 |
if (DoC[j + m] < VMAXINT32) |
1684 |
gpertea |
29 |
DoC[j + m]++; |
1685 |
|
|
} |
1686 |
|
|
//fall through this case to REF_SKIP is intentional |
1687 |
|
|
case REF_SKIP: |
1688 |
|
|
j += cigar[c].length; |
1689 |
|
|
break; |
1690 |
|
|
default: |
1691 |
|
|
break; |
1692 |
|
|
} |
1693 |
|
|
|
1694 |
|
|
} |
1695 |
|
|
} |
1696 |
|
|
|
1697 |
|
|
void print_hit(FILE* fout, |
1698 |
|
|
const char* read_name, |
1699 |
|
|
const BowtieHit& bh, |
1700 |
|
|
const char* ref_name, |
1701 |
|
|
const char* sequence, |
1702 |
|
|
const char* qualities, |
1703 |
|
|
bool from_bowtie) |
1704 |
|
|
{ |
1705 |
|
|
string seq; |
1706 |
|
|
string quals; |
1707 |
|
|
if (sequence) |
1708 |
|
|
{ |
1709 |
|
|
seq = sequence; |
1710 |
|
|
quals = qualities; |
1711 |
|
|
seq.resize(bh.read_len()); |
1712 |
|
|
quals.resize(bh.read_len()); |
1713 |
|
|
} |
1714 |
|
|
else |
1715 |
|
|
{ |
1716 |
|
|
seq = "*"; |
1717 |
|
|
} |
1718 |
|
|
|
1719 |
|
|
if (qualities) |
1720 |
|
|
{ |
1721 |
|
|
quals = qualities; |
1722 |
|
|
quals.resize(bh.read_len()); |
1723 |
|
|
} |
1724 |
|
|
else |
1725 |
|
|
{ |
1726 |
|
|
quals = "*"; |
1727 |
|
|
} |
1728 |
|
|
|
1729 |
|
|
uint32_t sam_flag = 0; |
1730 |
|
|
if (bh.antisense_align()) |
1731 |
|
|
{ |
1732 |
|
|
sam_flag |= 0x0010; // BAM_FREVERSE |
1733 |
|
|
if (sequence && !from_bowtie) // if it is from bowtie hit, it's already reversed. |
1734 |
|
|
{ |
1735 |
|
|
reverse_complement(seq); |
1736 |
|
|
reverse(quals.begin(), quals.end()); |
1737 |
|
|
} |
1738 |
|
|
} |
1739 |
|
|
|
1740 |
|
|
uint32_t sam_pos = bh.left() + 1; |
1741 |
|
|
uint32_t map_quality = 255; |
1742 |
|
|
char cigar[256]; |
1743 |
|
|
cigar[0] = 0; |
1744 |
|
|
string mate_ref_name = "*"; |
1745 |
|
|
uint32_t mate_pos = 0; |
1746 |
|
|
uint32_t insert_size = 0; |
1747 |
|
|
//string qualities = "*"; |
1748 |
|
|
|
1749 |
|
|
const vector<CigarOp>& bh_cigar = bh.cigar(); |
1750 |
|
|
|
1751 |
|
|
/* |
1752 |
|
|
* In addition to calculating the cigar string, |
1753 |
|
|
* we need to figure out how many in/dels are in the |
1754 |
|
|
* sequence, so that we can give the correct |
1755 |
|
|
* value for the NM tag |
1756 |
|
|
*/ |
1757 |
|
|
int indel_distance = 0; |
1758 |
|
|
for (size_t c = 0; c < bh_cigar.size(); ++c) |
1759 |
|
|
{ |
1760 |
|
|
char ibuf[64]; |
1761 |
|
|
sprintf(ibuf, "%d", bh_cigar[c].length); |
1762 |
|
|
switch(bh_cigar[c].opcode) |
1763 |
|
|
{ |
1764 |
|
|
case MATCH: |
1765 |
|
|
strcat(cigar, ibuf); |
1766 |
|
|
strcat(cigar, "M"); |
1767 |
|
|
break; |
1768 |
|
|
case INS: |
1769 |
|
|
strcat(cigar, ibuf); |
1770 |
|
|
strcat(cigar, "I"); |
1771 |
|
|
indel_distance += bh_cigar[c].length; |
1772 |
|
|
break; |
1773 |
|
|
case DEL: |
1774 |
|
|
strcat(cigar, ibuf); |
1775 |
|
|
strcat(cigar, "D"); |
1776 |
|
|
indel_distance += bh_cigar[c].length; |
1777 |
|
|
break; |
1778 |
|
|
case REF_SKIP: |
1779 |
|
|
strcat(cigar, ibuf); |
1780 |
|
|
strcat(cigar, "N"); |
1781 |
|
|
break; |
1782 |
|
|
default: |
1783 |
|
|
break; |
1784 |
|
|
} |
1785 |
|
|
} |
1786 |
|
|
|
1787 |
|
|
//string q = string(bh.read_len(), '!'); |
1788 |
|
|
//string s = string(bh.read_len(), 'N'); |
1789 |
|
|
|
1790 |
|
|
fprintf(fout, |
1791 |
|
|
"%s\t%d\t%s\t%d\t%d\t%s\t%s\t%d\t%d\t%s\t%s", |
1792 |
|
|
read_name, |
1793 |
|
|
sam_flag, |
1794 |
|
|
ref_name, |
1795 |
|
|
sam_pos, |
1796 |
|
|
map_quality, |
1797 |
|
|
cigar, |
1798 |
|
|
mate_ref_name.c_str(), |
1799 |
|
|
mate_pos, |
1800 |
|
|
insert_size, |
1801 |
|
|
seq.c_str(), |
1802 |
|
|
quals.c_str()); |
1803 |
|
|
|
1804 |
|
|
if (!sam_readgroup_id.empty()) |
1805 |
|
|
fprintf(fout, "\tRG:Z:%s", sam_readgroup_id.c_str()); |
1806 |
|
|
|
1807 |
|
|
fprintf(fout, "\tNM:i:%d", bh.edit_dist() + indel_distance); |
1808 |
|
|
|
1809 |
|
|
bool containsSplice = false; |
1810 |
|
|
for(vector<CigarOp>::const_iterator itr = bh.cigar().begin(); itr != bh.cigar().end(); ++itr){ |
1811 |
|
|
if(itr->opcode == REF_SKIP){ |
1812 |
|
|
containsSplice = true; |
1813 |
|
|
break; |
1814 |
|
|
} |
1815 |
|
|
} |
1816 |
|
|
|
1817 |
|
|
if (containsSplice) |
1818 |
|
|
fprintf(fout, "\tXS:A:%c", bh.antisense_splice() ? '-' : '+'); |
1819 |
|
|
|
1820 |
|
|
fprintf(fout, "\n"); |
1821 |
|
|
} |
1822 |
|
|
|
1823 |
|
|
void print_bamhit(GBamWriter& wbam, |
1824 |
|
|
const char* read_name, |
1825 |
|
|
const BowtieHit& bh, |
1826 |
|
|
const char* ref_name, |
1827 |
|
|
const char* sequence, |
1828 |
|
|
const char* qualities, |
1829 |
|
|
bool from_bowtie) |
1830 |
|
|
{ |
1831 |
|
|
string seq; |
1832 |
|
|
string quals; |
1833 |
|
|
if (sequence) { |
1834 |
|
|
seq = sequence; |
1835 |
|
|
quals = qualities; |
1836 |
|
|
seq.resize(bh.read_len()); |
1837 |
|
|
quals.resize(bh.read_len()); |
1838 |
|
|
} |
1839 |
|
|
else { |
1840 |
|
|
seq = "*"; |
1841 |
|
|
} |
1842 |
|
|
if (qualities) { |
1843 |
|
|
quals = qualities; |
1844 |
|
|
quals.resize(bh.read_len()); |
1845 |
|
|
} |
1846 |
|
|
else |
1847 |
|
|
{ |
1848 |
|
|
quals = "*"; |
1849 |
|
|
} |
1850 |
|
|
|
1851 |
|
|
uint32_t sam_flag = 0; |
1852 |
|
|
if (bh.antisense_align()) |
1853 |
|
|
{ |
1854 |
|
|
sam_flag |= 0x0010; // BAM_FREVERSE |
1855 |
|
|
if (sequence && !from_bowtie) // if it is from bowtie hit, it's already reversed. |
1856 |
|
|
{ |
1857 |
|
|
reverse_complement(seq); |
1858 |
|
|
reverse(quals.begin(), quals.end()); |
1859 |
|
|
} |
1860 |
|
|
} |
1861 |
|
|
|
1862 |
|
|
uint32_t sam_pos = bh.left() + 1; |
1863 |
|
|
uint32_t map_quality = 255; |
1864 |
|
|
char cigar[256]; |
1865 |
|
|
cigar[0] = 0; |
1866 |
|
|
string mate_ref_name = "*"; |
1867 |
|
|
uint32_t mate_pos = 0; |
1868 |
|
|
uint32_t insert_size = 0; |
1869 |
|
|
//string qualities = "*"; |
1870 |
|
|
|
1871 |
|
|
const vector<CigarOp>& bh_cigar = bh.cigar(); |
1872 |
|
|
/* |
1873 |
|
|
* In addition to calculating the cigar string, |
1874 |
|
|
* we need to figure out how many in/dels are in the |
1875 |
|
|
* sequence, so that we can give the correct |
1876 |
|
|
* value for the NM tag |
1877 |
|
|
*/ |
1878 |
|
|
int indel_distance = 0; |
1879 |
|
|
for (size_t c = 0; c < bh_cigar.size(); ++c) |
1880 |
|
|
{ |
1881 |
|
|
char ibuf[64]; |
1882 |
|
|
sprintf(ibuf, "%d", bh_cigar[c].length); |
1883 |
|
|
switch(bh_cigar[c].opcode) |
1884 |
|
|
{ |
1885 |
|
|
case MATCH: |
1886 |
|
|
strcat(cigar, ibuf); |
1887 |
|
|
strcat(cigar, "M"); |
1888 |
|
|
break; |
1889 |
|
|
case INS: |
1890 |
|
|
strcat(cigar, ibuf); |
1891 |
|
|
strcat(cigar, "I"); |
1892 |
|
|
indel_distance += bh_cigar[c].length; |
1893 |
|
|
break; |
1894 |
|
|
case DEL: |
1895 |
|
|
strcat(cigar, ibuf); |
1896 |
|
|
strcat(cigar, "D"); |
1897 |
|
|
indel_distance += bh_cigar[c].length; |
1898 |
|
|
break; |
1899 |
|
|
case REF_SKIP: |
1900 |
|
|
strcat(cigar, ibuf); |
1901 |
|
|
strcat(cigar, "N"); |
1902 |
|
|
break; |
1903 |
|
|
default: |
1904 |
|
|
break; |
1905 |
|
|
} |
1906 |
|
|
} |
1907 |
|
|
|
1908 |
|
|
//string q = string(bh.read_len(), '!'); |
1909 |
|
|
//string s = string(bh.read_len(), 'N'); |
1910 |
|
|
/* |
1911 |
|
|
fprintf(fout, |
1912 |
|
|
"%s\t%d\t%s\t%d\t%d\t%s\t%s\t%d\t%d\t%s\t%s", |
1913 |
|
|
read_name, |
1914 |
|
|
sam_flag, |
1915 |
|
|
ref_name, |
1916 |
|
|
sam_pos, |
1917 |
|
|
map_quality, |
1918 |
|
|
cigar, |
1919 |
|
|
mate_ref_name.c_str(), |
1920 |
|
|
mate_pos, |
1921 |
|
|
insert_size, |
1922 |
|
|
seq.c_str(), |
1923 |
|
|
quals.c_str()); |
1924 |
|
|
|
1925 |
|
|
fprintf(fout, "\tNM:i:%d", bh.edit_dist() + indel_distance); |
1926 |
|
|
*/ |
1927 |
|
|
vector<string> auxdata; |
1928 |
|
|
|
1929 |
|
|
if (!sam_readgroup_id.empty()) |
1930 |
|
|
{ |
1931 |
|
|
string nm("RG:Z:"); |
1932 |
|
|
nm += sam_readgroup_id; |
1933 |
|
|
auxdata.push_back(nm); |
1934 |
|
|
} |
1935 |
|
|
|
1936 |
|
|
string nm("NM:i:"); |
1937 |
|
|
str_appendInt(nm, bh.edit_dist() + indel_distance); |
1938 |
|
|
auxdata.push_back(nm); |
1939 |
|
|
bool containsSplice = false; |
1940 |
|
|
for(vector<CigarOp>::const_iterator itr = bh.cigar().begin(); itr != bh.cigar().end(); ++itr) |
1941 |
|
|
if(itr->opcode == REF_SKIP){ |
1942 |
|
|
containsSplice = true; |
1943 |
|
|
break; |
1944 |
|
|
} |
1945 |
|
|
|
1946 |
|
|
if (containsSplice) { |
1947 |
|
|
//fprintf(fout, "\tXS:A:%c", bh.antisense_splice() ? '-' : '+'); |
1948 |
|
|
nm="XS:A:"; |
1949 |
|
|
nm+=(char)(bh.antisense_splice() ? '-' : '+'); |
1950 |
|
|
auxdata.push_back(nm); |
1951 |
|
|
} |
1952 |
|
|
|
1953 |
|
|
GBamRecord *brec = wbam.new_record(read_name, sam_flag, ref_name, sam_pos, map_quality, |
1954 |
|
|
cigar, mate_ref_name.c_str(), mate_pos, |
1955 |
|
|
insert_size, seq.c_str(), quals.c_str(), &auxdata); |
1956 |
|
|
|
1957 |
|
|
|
1958 |
|
|
|
1959 |
|
|
wbam.write(brec); |
1960 |
|
|
delete brec; |
1961 |
|
|
} |
1962 |
|
|
|
1963 |
|
|
/** |
1964 |
|
|
* Print a vector of cigar operations to a file. |
1965 |
|
|
* @param bh_cigar A vector of CigarOps |
1966 |
|
|
* @return a string representation of the cigar string |
1967 |
|
|
*/ |
1968 |
|
|
std::string print_cigar(vector<CigarOp>& bh_cigar){ |
1969 |
|
|
char cigar[256]; |
1970 |
|
|
cigar[0] = 0; |
1971 |
|
|
for (size_t c = 0; c < bh_cigar.size(); ++c) |
1972 |
|
|
{ |
1973 |
|
|
char ibuf[64]; |
1974 |
|
|
sprintf(ibuf, "%d", bh_cigar[c].length); |
1975 |
|
|
switch(bh_cigar[c].opcode) |
1976 |
|
|
{ |
1977 |
|
|
case MATCH: |
1978 |
|
|
strcat(cigar, ibuf); |
1979 |
|
|
strcat(cigar, "M"); |
1980 |
|
|
break; |
1981 |
|
|
case INS: |
1982 |
|
|
strcat(cigar, ibuf); |
1983 |
|
|
strcat(cigar, "I"); |
1984 |
|
|
break; |
1985 |
|
|
case DEL: |
1986 |
|
|
strcat(cigar, ibuf); |
1987 |
|
|
strcat(cigar, "D"); |
1988 |
|
|
break; |
1989 |
|
|
case REF_SKIP: |
1990 |
|
|
strcat(cigar, ibuf); |
1991 |
|
|
strcat(cigar, "N"); |
1992 |
|
|
break; |
1993 |
|
|
default: |
1994 |
|
|
break; |
1995 |
|
|
} |
1996 |
|
|
} |
1997 |
|
|
string result(cigar); |
1998 |
|
|
return result; |
1999 |
|
|
} |
2000 |
|
|
|
2001 |
|
|
bool BowtieHit::check_editdist_consistency(const RefSequenceTable& rt) |
2002 |
|
|
{ |
2003 |
|
|
RefSequenceTable::Sequence* ref_str = rt.get_seq(_ref_id); |
2004 |
|
|
if (!ref_str) |
2005 |
|
|
return false; |
2006 |
|
|
|
2007 |
|
|
const seqan::Dna5String ref_seq = seqan::infix(*ref_str, _left, right()); |
2008 |
|
|
|
2009 |
|
|
size_t pos_seq = 0; |
2010 |
|
|
size_t pos_ref = 0; |
2011 |
|
|
size_t mismatch = 0; |
2012 |
|
|
for (size_t i = 0; i < _cigar.size(); ++i) |
2013 |
|
|
{ |
2014 |
|
|
CigarOp cigar = _cigar[i]; |
2015 |
|
|
switch(cigar.opcode) |
2016 |
|
|
{ |
2017 |
|
|
case MATCH: |
2018 |
|
|
{ |
2019 |
|
|
for (size_t j = 0; j < cigar.length; ++j) |
2020 |
|
|
{ |
2021 |
|
|
if (_seq[pos_seq++] != ref_seq[pos_ref++]) |
2022 |
|
|
++mismatch; |
2023 |
|
|
} |
2024 |
|
|
} |
2025 |
|
|
break; |
2026 |
|
|
case INS: |
2027 |
|
|
{ |
2028 |
|
|
pos_seq += cigar.length; |
2029 |
|
|
} |
2030 |
|
|
break; |
2031 |
|
|
|
2032 |
|
|
case DEL: |
2033 |
|
|
case REF_SKIP: |
2034 |
|
|
{ |
2035 |
|
|
pos_ref += cigar.length; |
2036 |
|
|
} |
2037 |
|
|
break; |
2038 |
|
|
|
2039 |
|
|
default: |
2040 |
|
|
break; |
2041 |
|
|
} |
2042 |
|
|
} |
2043 |
|
|
|
2044 |
|
|
return mismatch == _edit_dist; |
2045 |
|
|
} |