ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/tophat_cpp/segment_juncs.cpp
Revision: 29
Committed: Tue Aug 2 21:24:54 2011 UTC (13 years, 1 month ago) by gpertea
File size: 119886 byte(s)
Log Message: 
adding tophat source work
Line File contents
1 #ifdef HAVE_CONFIG_H
2 #include <config.h>
3 #endif
4
5 /*
6 * segment_juncs.cpp
7 * TopHat
8 *
9 * Created by Cole Trapnell on 2/5/09.
10 * Copyright 2009 Cole Trapnell. All rights reserved.
11 *
12 */
13
14 #include <cassert>
15 #include <cstdio>
16 #include <vector>
17 #include <string>
18 #include <map>
19 #include <algorithm>
20 #include <set>
21 //#include <stdexcept>
22 #include <iostream>
23 #include <fstream>
24 #include <sstream>
25 #include <cstring>
26 #include <bitset>
27 #include <seqan/sequence.h>
28 #include <seqan/find.h>
29 #include <seqan/file.h>
30 #include <seqan/modifier.h>
31 #include <getopt.h>
32
33 #include "common.h"
34 #include "bwt_map.h"
35 #include "tokenize.h"
36 #include "segments.h"
37 #include "reads.h"
38 #include "junctions.h"
39 #include "insertions.h"
40 #include "deletions.h"
41
42 using namespace seqan;
43 using namespace std;
44 using namespace __gnu_cxx;
45
46 // daehwan
47 bool bDebug = false;
48
49 void print_usage()
50 {
51 fprintf(stderr, "Usage: segment_juncs <ref.fa> <segment.juncs> <segment.insertions> <segment.deletions> <left_reads.fq> <left_seg1.bwtout,...,segN.bwtout> [right_reads.fq right_seg1.bwtout,...,right_segN.bwtout]\n");
52 }
53
54 // This is the maximum number of bowtie mismatches allower per segment hit
55 static const int num_bowtie_mismatches = 2;
56
57 static const int max_cov_juncs = 5000000;
58 static const int max_seg_juncs = 10000000;
59 int max_microexon_stretch = 2000;
60 int butterfly_overhang = 6;
61
62 void get_seqs(istream& ref_stream,
63 RefSequenceTable& rt,
64 bool keep_seqs = true,
65 bool strip_slash = false)
66 {
67 while(ref_stream.good() &&
68 !ref_stream.eof())
69 {
70 RefSequenceTable::Sequence* ref_str = new RefSequenceTable::Sequence();
71 string name;
72 readMeta(ref_stream, name, Fasta());
73 string::size_type space_pos = name.find_first_of(" \t\r");
74 if (space_pos != string::npos)
75 {
76 name.resize(space_pos);
77 }
78 fprintf(stderr, "\tLoading %s...", name.c_str());
79 seqan::read(ref_stream, *ref_str, Fasta());
80 fprintf(stderr, "done\n");
81 rt.get_id(name, keep_seqs ? ref_str : NULL, 0);
82 if (!keep_seqs)
83 delete ref_str;
84 }
85 }
86
87
88 RefSeg seg_from_bowtie_hit(const BowtieHit& T)
89 {
90 RefSeg r_seg(T.ref_id(), POINT_DIR_DONTCARE, T.antisense_align(), READ_DONTCARE, 0, 0);
91
92 if (T.antisense_align())
93 {
94 r_seg.left = max(0, T.right() - 2);
95 r_seg.right = T.right() + (T.right() - T.left() + 1); // num allowed bowtie mismatches
96 r_seg.points_where = POINT_DIR_RIGHT;
97 }
98 else
99 {
100 r_seg.left = max(0, T.left() - (T.right() - T.left() + 1));
101 r_seg.right = T.left() + 2; // num allowed bowtie mismatches
102 r_seg.points_where = POINT_DIR_LEFT;
103 }
104
105 return r_seg;
106 }
107
108 pair<RefSeg, RefSeg> segs_from_bowtie_hits(const BowtieHit& T,
109 const BowtieHit& H)
110 {
111 pair<RefSeg, RefSeg> seg_pair;
112 if (H.antisense_align() == false &&
113 abs((H.right() + 1) - T.left()) < (int)max_segment_intron_length)
114 {
115 RefSeg left_seg(H.ref_id(), POINT_DIR_RIGHT, H.antisense_align(), READ_DONTCARE, 0, 0);
116 left_seg.left = max(0, H.right() - 2);
117 left_seg.right = H.right() + (H.right() - H.left() + 1); // num allowed bowtie mismatches
118
119 RefSeg right_seg(T.ref_id(), POINT_DIR_LEFT, T.antisense_align(), READ_DONTCARE, 0, 0);
120 right_seg.left = max(0, T.left() - (T.right() - T.left() + 1));
121 right_seg.right = T.left() + 2; // num allowed bowtie mismatches
122
123 seg_pair = make_pair(left_seg, right_seg);
124 }
125 else if (H.antisense_align() == true &&
126 abs((T.right() + 1) - H.left()) < (int)max_segment_intron_length)
127 {
128 RefSeg left_seg(T.ref_id(), POINT_DIR_RIGHT, T.antisense_align(), READ_DONTCARE, 0, 0);
129 left_seg.left = max(0, T.right() - 2);
130 left_seg.right = T.right() + (T.right() - T.left() + 1); // num allowed bowtie mismatches
131
132 RefSeg right_seg(H.ref_id(), POINT_DIR_LEFT, H.antisense_align(), READ_DONTCARE, 0, 0);
133 right_seg.left = max(0, H.left() - (H.right() - H.left() + 1));
134 right_seg.right = H.left() + 2; // num allowed bowtie mismatches
135
136 seg_pair = make_pair(left_seg, right_seg);
137 }
138 return seg_pair;
139 }
140
141 //static const size_t half_splice_mer_len = 6;
142 //static const size_t splice_mer_len = 2 * half_splice_mer_len;
143
144 struct MerExtension
145 {
146 static const int MAX_EXTENSION_BP = 14;
147 uint32_t left_dna_str : 28; // up to 14bp encoded in 2-bits-per-base
148 uint8_t left_ext_len : 4; // how many bases in this extension
149
150 uint32_t right_dna_str : 28; // up to 14bp encoded in 2-bits-per-base
151 uint8_t right_ext_len : 4; // how many bases in this extension
152
153 MerExtension() : left_dna_str(0), left_ext_len(0), right_dna_str(0), right_ext_len(0) {}
154
155 bool operator<(const MerExtension& rhs) const
156 {
157 if (left_dna_str != rhs.left_dna_str)
158 return left_dna_str < rhs.left_dna_str;
159 if (left_ext_len != rhs.left_ext_len)
160 return left_ext_len < rhs.left_ext_len;
161 if (right_dna_str != rhs.right_dna_str)
162 return right_dna_str < rhs.right_dna_str;
163 if (right_ext_len != rhs.right_ext_len)
164 return right_ext_len < rhs.right_ext_len;
165 return false;
166 }
167
168 bool operator==(const MerExtension& rhs) const
169 {
170 bool eq = left_dna_str == rhs.left_dna_str &&
171 right_dna_str == rhs.right_dna_str &&
172 left_ext_len == rhs.left_ext_len &&
173 right_ext_len == rhs.right_ext_len;
174
175 return eq;
176 }
177 };
178
179
180 /// For converting from ASCII to the Dna5 code where A=0, C=1, G=2,
181 /// T=3, N=4
182 uint8_t charToDna5[] = {
183 /* 0 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
184 /* 16 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
185 /* 32 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
186 /* 48 */ 0, 1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
187 /* 64 */ 0, 0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 4, 0,
188 /* A C G N */
189 /* 80 */ 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
190 /* T */
191 /* 96 */ 0, 0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 4, 0,
192 /* a c g n */
193 /* 112 */ 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
194 /* t */
195 /* 128 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
196 /* 144 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
197 /* 160 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
198 /* 176 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
199 /* 192 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
200 /* 208 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
201 /* 224 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
202 /* 240 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
203 };
204
205
206 typedef vector<vector<MerExtension> > MerExtensionTable;
207 typedef vector<uint32_t> MerExtensionCounts;
208
209 MerExtensionTable extensions;
210 MerExtensionCounts extension_counts;
211
212 uint64_t dna5str_to_idx(const string& str)
213 {
214 uint64_t idx = 0;
215
216 for (size_t i = 0; i < str.length(); ++i)
217 {
218 idx <<=2;
219 char c = toupper(str[i]);
220 idx |= (0x3 & charToDna5[(size_t)c]);
221 }
222 return idx;
223 }
224
225 uint64_t colorstr_to_idx(const string& str)
226 {
227 uint64_t idx = 0;
228
229 for (size_t i = 0; i < str.length(); ++i)
230 {
231 idx <<=2;
232 char c = str[i];
233 idx |= (0x3 & charToDna5[(size_t)c]);
234 }
235 return idx;
236 }
237
238 void store_read_extensions(MerExtensionTable& ext_table,
239 int seq_key_len,
240 int min_ext_len,
241 const string& seq,
242 bool use_precount_table)
243 {
244 // h is will hold the 2-bit-per-base representation of the k-mer seeds for
245 // this read.
246
247 uint64_t seed = 0;
248 bitset<256> left = 0;
249 bitset<256> right = 0;
250 const char* p = seq.c_str();
251
252 unsigned int seq_len = (int)seq.length();
253 const char* seq_end = p + seq_len;
254
255 // Build the first seed
256 while (p < seq.c_str() + (2 * seq_key_len))
257 {
258 seed <<= 2;
259 seed |= (0x3 & charToDna5[(size_t)*p]);
260 ++p;
261 }
262
263 // Build the rest of them with a sliding window, adding successive bases
264 // to the "right-remainder" word.
265 while (p < seq_end)
266 {
267 right <<= 2;
268 right |= (0x3 & charToDna5[(size_t)*p]);
269 ++p;
270 }
271
272 // This loop will construct successive seed, along with 32-bit words
273 // containing the left and right remainders for each seed
274 uint32_t i = 0;
275 size_t new_hits = 0;
276 do
277 {
278 // How many base pairs exist in the right remainder beyond what we have
279 // space for ?
280 int extra_right_bp = ((int)seq.length() -
281 (i + 2 * seq_key_len)) - MerExtension::MAX_EXTENSION_BP;
282
283 uint32_t hit_right = 0;
284 if (extra_right_bp > 0)
285 {
286 //bitset<32> tmp_hit_right = (right >> (extra_right_bp << 1));
287 hit_right = (uint32_t)(right >> (extra_right_bp << 1)).to_ulong();
288 }
289 else
290 {
291 hit_right = (uint32_t)right.to_ulong();
292 }
293
294 uint32_t hit_left = (uint32_t)((left << (256 - 32)) >> (256 - 32)).to_ulong();
295
296 //size_t prev_cap = (*mer_table)[seed].capacity();
297 //(*mer_table)[seed].push_back(ReadHit(hit_left, hit_right,i, read_num, reverse_complement));
298 //cap_increase += ((*mer_table)[seed].capacity() - prev_cap) * sizeof (ReadHit);
299
300 MerExtension ext;
301
302 ext.right_dna_str = hit_right;
303 ext.right_ext_len = min(seq_len - (2 * seq_key_len) - i,
304 (unsigned int)MerExtension::MAX_EXTENSION_BP);
305
306 ext.left_dna_str = hit_left;
307 ext.left_ext_len = min(i, (unsigned int)MerExtension::MAX_EXTENSION_BP);
308 if (use_precount_table)
309 {
310 int curr_seed = --extension_counts[seed];
311 if (curr_seed < 0 || curr_seed > (int)ext_table[seed].size())
312 {
313 fprintf(stderr, "Error: curr_seed is %d, max is %lu\n", curr_seed, (long unsigned int)ext_table[seed].size());
314 }
315
316 ext_table[seed][curr_seed] = ext;
317 }
318 else
319 {
320 ext_table[seed].push_back(ext);
321 }
322 new_hits++;
323
324 // Take the leftmost base of the seed and stick it into bp
325 uint64_t bp = seed & (0x3uLL << ((seq_key_len << 2) - 2));
326
327 // Move that base down to the least significant bits of bp
328 bp >>= ((seq_key_len << 2) - 2);
329
330 // And tack it onto the left remainder of the read
331 left <<= 2;
332 left |= bp;
333
334 // Now take the leftmost base of the right remainder and stick it into
335 // the rightmost position of the seed
336 uint32_t right_len = seq_len - (i + seq_key_len * 2);
337 //bp = right & (0x3uLL << ((right_len - 1) << 1));
338
339 seed <<= 2;
340 //cout << right << endl;
341 bitset<256> tmp_right = (right >> ((right_len - 1) << 1));
342 //cout <<tmp_right << endl;
343 seed |= tmp_right.to_ulong();
344 seed &= ~(0xFFFFFFFFFFFFFFFFuLL << (seq_key_len << 2));
345
346
347 //Now remove that leftmost base of the right remainder
348 //right &= ~(0xFFFFFFFFFFFFFFFFuLL << (right_len - 1 << 1));
349 if (right_len)
350 {
351 right.set(((right_len - 1) << 1), 0);
352 right.set(((right_len - 1) << 1) + 1, 0);
353 }
354 ++i;
355
356 }while(i <= (size_t)(seq_end - seq.c_str()) - (2 * seq_key_len));
357 return;
358 }
359
360 void count_read_extensions(MerExtensionCounts& ext_counts,
361 int seq_key_len,
362 int min_ext_len,
363 const string& seq)
364 {
365 // h is will hold the 2-bit-per-base representation of the k-mer seeds for
366 // this read.
367
368 uint64_t seed = 0;
369 bitset<256> left = 0;
370 bitset<256> right = 0;
371 const char* p = seq.c_str();
372
373 unsigned int seq_len = (int)seq.length();
374 const char* seq_end = p + seq_len;
375
376 // Build the first seed
377 while (p < seq.c_str() + (2 * seq_key_len))
378 {
379 seed <<= 2;
380 seed |= (0x3 & charToDna5[(size_t)*p]);
381 ++p;
382 }
383
384 // Build the rest of them with a sliding window, adding successive bases
385 // to the "right-remainder" word.
386 while (p < seq_end)
387 {
388 right <<= 2;
389 right |= (0x3 & charToDna5[(size_t)*p]);
390 ++p;
391 }
392
393 // This loop will construct successive seed, along with 32-bit words
394 // containing the left and right remainders for each seed
395 uint32_t i = 0;
396 size_t new_hits = 0;
397 do
398 {
399 ext_counts[seed]++;
400
401 new_hits++;
402
403 // Take the leftmost base of the seed and stick it into bp
404 uint64_t bp = seed & (0x3uLL << ((seq_key_len << 2) - 2));
405
406 // Move that base down to the least significant bits of bp
407 bp >>= ((seq_key_len << 2) - 2);
408
409 // And tack it onto the left remainder of the read
410 left <<= 2;
411 left |= bp;
412
413 // Now take the leftmost base of the right remainder and stick it into
414 // the rightmost position of the seed
415 uint32_t right_len = seq_len - (i + seq_key_len * 2);
416 //bp = right & (0x3uLL << ((right_len - 1) << 1));
417
418 seed <<= 2;
419 //cout << right << endl;
420 bitset<256> tmp_right = (right >> ((right_len - 1) << 1));
421 //cout <<tmp_right << endl;
422 seed |= tmp_right.to_ulong();
423 seed &= ~(0xFFFFFFFFFFFFFFFFuLL << (seq_key_len << 2));
424
425
426 //Now remove that leftmost base of the right remainder
427 //right &= ~(0xFFFFFFFFFFFFFFFFuLL << (right_len - 1 << 1));
428 if (right_len)
429 {
430 right.set(((right_len - 1) << 1), 0);
431 right.set(((right_len - 1) << 1) + 1, 0);
432 }
433 ++i;
434
435 }while(i <= (size_t)(seq_end - seq.c_str()) - (2 * seq_key_len));
436 return;
437 }
438
439 //void count_read_mers(FILE* reads_file, size_t half_splice_mer_len)
440 void count_read_mers(FZPipe& reads_file, size_t half_splice_mer_len)
441 {
442 Read read;
443 size_t splice_mer_len = 2 * half_splice_mer_len;
444 size_t mer_table_size = 1 << ((splice_mer_len)<<1);
445 extension_counts.resize(mer_table_size);
446 FLineReader fr(reads_file);
447 //while(!feof(reads_file))
448 while (!fr.isEof())
449 {
450 read.clear();
451
452 // Get the next read from the file
453 if (!next_fasta_record(fr, read.name, read.seq, reads_format))
454 break;
455 if (reads_format == FASTQ)
456 {
457 if (!next_fastq_record(fr, read.seq, read.alt_name, read.qual, reads_format))
458 break;
459 }
460
461 if (color)
462 // erase the primer and the adjacent color
463 read.seq.erase(0, 2);
464
465 if (read.seq.size() > 32)
466 read.seq.resize(32);
467
468 count_read_extensions(extension_counts,
469 half_splice_mer_len,
470 half_splice_mer_len,
471 read.seq);
472 }
473
474 //rewind(reads_file);
475 reads_file.rewind();
476 }
477
478 void compact_extension_table()
479 {
480 for (size_t i = 0; i < extensions.size(); ++i)
481 {
482 vector<MerExtension>& exts = extensions[i];
483 sort(exts.begin(), exts.end());
484 vector<MerExtension>::iterator new_end = unique(exts.begin(), exts.end());
485 exts.erase(new_end, exts.end());
486 vector<MerExtension>(exts).swap(exts);
487 }
488 }
489
490 void prune_extension_table(uint8_t max_extension_bp)
491 {
492 uint32_t mask = ~(0xFFFFFFFFuLL << (max_extension_bp << 1));
493
494 for (size_t i = 0; i < extensions.size(); ++i)
495 {
496 vector<MerExtension>& exts = extensions[i];
497 for (size_t j = 0; j < exts.size(); ++j)
498 {
499 MerExtension& ex = exts[j];
500 if (ex.left_ext_len > max_extension_bp)
501 {
502 ex.left_ext_len = max_extension_bp;
503 ex.left_dna_str &= mask;
504 }
505
506 if (ex.right_ext_len > max_extension_bp)
507 {
508 ex.right_dna_str >>= ((ex.right_ext_len - max_extension_bp) << 1);
509 ex.right_ext_len = max_extension_bp;
510 }
511 }
512 }
513 }
514
515 //void store_read_mers(FILE* reads_file, size_t half_splice_mer_len)
516 void store_read_mers(FZPipe& reads_file, size_t half_splice_mer_len)
517 {
518 Read read;
519 size_t splice_mer_len = 2 * half_splice_mer_len;
520
521 size_t mer_table_size = 1 << ((splice_mer_len)<<1);
522 extensions.resize(mer_table_size);
523
524 size_t num_indexed_reads = 0;
525 FLineReader fr(reads_file);
526 //while(!feof(reads_file))
527 while(!fr.isEof())
528 {
529 read.clear();
530
531 // Get the next read from the file
532 if (!next_fasta_record(fr, read.name, read.seq, reads_format))
533 break;
534 if (reads_format == FASTQ)
535 {
536 if (!next_fastq_record(fr, read.seq, read.alt_name, read.qual, reads_format))
537 break;
538 }
539
540 if (color)
541 // erase the primer and the adjacent color
542 read.seq.erase(0, 2);
543
544 if (read.seq.size() > 32)
545 read.seq.resize(32);
546
547 store_read_extensions(extensions,
548 half_splice_mer_len,
549 half_splice_mer_len,
550 read.seq,
551 true);
552
553 // Do NOT index the reverse of the reads
554
555 num_indexed_reads++;
556 if (num_indexed_reads % 1000000 == 0)
557 {
558 //fprintf(stderr, "Indexed %lu reads, compacting extension table\n", num_indexed_reads);
559 //compact_extension_table();
560 }
561 }
562
563 //fprintf(stderr, "Indexed %lu reads, compacting extension table\n", num_indexed_reads)
564
565 uint64_t num_extensions = 0;
566 for (size_t i = 0; i < extensions.size(); ++i)
567 {
568 num_extensions += extensions[i].size();
569 }
570 fprintf (stderr, "Total extensions: %lu\n", (long unsigned int)num_extensions);
571 //rewind(reads_file);
572 reads_file.rewind();
573 }
574
575 //void index_read_mers(vector<FILE*> reads_files,
576 void index_read_mers(vector<FZPipe>& reads_files,
577 size_t half_splice_mer_len)
578 {
579 extensions.clear();
580 for (size_t i = 0; i < reads_files.size(); ++i)
581 {
582 count_read_mers(reads_files[i], half_splice_mer_len);
583 }
584
585 extensions.resize(extension_counts.size());
586
587 for (size_t i = 0; i < extension_counts.size(); ++i)
588 {
589 extensions[i].resize(extension_counts[i]);
590 }
591
592 for (size_t i = 0; i < reads_files.size(); ++i)
593 {
594 store_read_mers(reads_files[i], half_splice_mer_len);
595 }
596
597 compact_extension_table();
598 }
599
600 /** Returns the number of characters in strings w1 and w2 that match,
601 * starting from right to left
602 */
603 int get_matching_chars(uint32_t w1, uint32_t w2)
604 {
605 //find the least significant mismatching bit between w1 and w2
606 int mismatch_bit = ffs(w1 ^ w2);
607
608 // If there is no mismatching bit, the words are equal
609 if (!mismatch_bit)
610 return -1;
611
612 // Given the mismatching bit, determine where the mismatching base is
613 mismatch_bit -= 1;
614 mismatch_bit -= ((mismatch_bit) & 1);
615 mismatch_bit >>= 1;
616
617 // Return the number of matching characters.
618 return mismatch_bit;
619 }
620
621 /**
622 * Computes the Hamming distance between two 16bp words, up to a specified
623 * maximum number of mismatches.
624 */
625 uint32_t mismatching_bases(uint32_t w1_word,
626 uint32_t w2_word,
627 int len,
628 uint32_t max_mis)
629 {
630 uint32_t diffs = 0;
631
632 int shift = 0;
633 int L = len;
634 uint32_t misses = 0;
635
636 // While we haven't yet exceeded the maximum allowable mismatches,
637 // and there are still unaligned bases, keep shift-anding
638 while (shift < len && misses <= max_mis)
639 {
640 int match_chars = 0;
641
642 // Get the number of characters matching on the right sides of
643 // both words
644 match_chars = get_matching_chars(w1_word, w2_word);
645
646 // If they are equal for this shift, we are done,
647 // the loop will stop at the next iteration
648 if (match_chars == -1 || match_chars >= len)
649 {
650 match_chars = len;
651 shift = len;
652 }
653 else
654 {
655 // If there is a mismatch in the remaining words
656 // decide how much to shift by and try again
657 match_chars = min(len, match_chars);
658 int shift_chars = (match_chars + 1);
659
660 L -= shift_chars;
661
662 int shift_bits = shift_chars << 1;
663
664 // Shift right past the matching part and the first mismatching base
665 w1_word >>= (shift_bits);
666 w2_word >>= (shift_bits);
667
668 shift += shift_chars;
669 diffs++;
670 misses++;
671 }
672 }
673
674 return diffs;
675 }
676
677 uint64_t rc_dna_str(uint64_t dna_str)
678 {
679 dna_str = ~dna_str;
680 uint64_t rc = 0;
681 for (int i = 0; i < 32; i++)
682 {
683 rc <<= 2;
684 rc |= (dna_str & 0x3);
685 dna_str >>= 2;
686 }
687 return rc;
688 }
689
690 uint64_t rc_color_str(uint64_t color_str)
691 {
692 uint64_t rc = 0;
693 for (int i = 0; i < 32; ++i)
694 {
695 rc <<= 2;
696 rc |= (color_str & 0x3);
697 color_str >>= 2;
698 }
699 return rc;
700 }
701
702 struct DnaSpliceStrings
703 {
704 DnaSpliceStrings(uint64_t f, uint64_t r) : fwd_string(f), rev_string(r), first_in_string('N'), last_in_string('N') {}
705 uint64_t fwd_string;
706 uint64_t rev_string;
707
708 // for color-space purposes
709 char first_in_string;
710 char last_in_string;
711
712 bool operator<(const DnaSpliceStrings& rhs) const
713 {
714 if (fwd_string != rhs.fwd_string)
715 return fwd_string < rhs.fwd_string;
716 if (rev_string != rhs.rev_string)
717 return rev_string < rhs.rev_string;
718 return false;
719 }
720
721 bool operator==(const DnaSpliceStrings& rhs) const
722 {
723 return fwd_string == rhs.fwd_string && rev_string == rhs.rev_string;
724 }
725 };
726
727 struct IntronMotifs
728 {
729 IntronMotifs(uint32_t rid) : ref_id(rid) {}
730 uint32_t ref_id;
731
732 vector<pair<size_t, DnaSpliceStrings> > fwd_donors;
733 vector<pair<size_t, DnaSpliceStrings> > fwd_acceptors;
734 vector<pair<size_t, DnaSpliceStrings> > rev_donors;
735 vector<pair<size_t, DnaSpliceStrings> > rev_acceptors;
736
737 void unique(vector<pair<size_t, DnaSpliceStrings> >& f)
738 {
739 sort(f.begin(), f.end());
740 vector<pair<size_t, DnaSpliceStrings> >::iterator i = std::unique(f.begin(), f.end());
741 f.erase(i, f.end());
742 }
743
744 void unique()
745 {
746 unique(fwd_donors);
747 unique(fwd_acceptors);
748 unique(rev_donors);
749 unique(rev_acceptors);
750 }
751
752 void attach_mers(RefSequenceTable::Sequence& ref_str)
753 {
754 attach_upstream_mers(ref_str, fwd_donors);
755 attach_upstream_mers(ref_str, rev_acceptors);
756
757 attach_downstream_mers(ref_str, rev_donors);
758 attach_downstream_mers(ref_str, fwd_acceptors);
759 }
760
761 void attach_upstream_mers(RefSequenceTable::Sequence& ref_str,
762 vector<pair<size_t, DnaSpliceStrings> >& dinucs)
763 {
764 for (size_t i = 0; i < dinucs.size(); ++i)
765 {
766 size_t pos = dinucs[i].first;
767 int half_splice_mer_len = 32;
768
769 if (color)
770 {
771 if (pos <= (size_t)half_splice_mer_len+1 || pos >= length(ref_str))
772 continue;
773
774 Dna5String seg_str = seqan::infixWithLength(ref_str,
775 pos - half_splice_mer_len - 1,
776 half_splice_mer_len + 1);
777 stringstream ss(stringstream::in | stringstream::out);
778 string s;
779 ss << seg_str;
780 ss >> s;
781
782 string col_seg_str = convert_bp_to_color(s, true);
783 uint64_t idx = colorstr_to_idx(col_seg_str);
784
785 dinucs[i].second.fwd_string = idx;
786 dinucs[i].second.rev_string = rc_color_str(idx);
787 dinucs[i].second.first_in_string = s[1];
788 dinucs[i].second.last_in_string = s[half_splice_mer_len];
789 }
790 else
791 {
792 if (pos <= (size_t)half_splice_mer_len || pos >= length(ref_str))
793 continue;
794
795 Dna5String seg_str = seqan::infixWithLength(ref_str,
796 pos - half_splice_mer_len,
797 half_splice_mer_len);
798
799 stringstream ss(stringstream::in | stringstream::out);
800 string s;
801 ss << seg_str;
802 ss >> s;
803 uint64_t idx = dna5str_to_idx(s);
804 dinucs[i].second.fwd_string = idx;
805 dinucs[i].second.rev_string = rc_dna_str(idx);
806 }
807 }
808 }
809
810
811 void attach_downstream_mers(RefSequenceTable::Sequence& ref_str,
812 vector<pair<size_t, DnaSpliceStrings> >& dinucs)
813 {
814 for (size_t i = 0; i < dinucs.size(); ++i)
815 {
816 size_t pos = dinucs[i].first;
817
818 int half_splice_mer_len = 32;
819
820 if (pos + 2 + half_splice_mer_len >= length(ref_str))
821 continue;
822
823 if (color)
824 {
825 Dna5String seg_str = seqan::infixWithLength(ref_str,
826 pos + 2 - 1,
827 half_splice_mer_len + 1);
828 stringstream ss(stringstream::in | stringstream::out);
829 string s;
830 ss << seg_str;
831 ss >> s;
832
833 string col_seg_str = convert_bp_to_color(s, true);
834 uint64_t idx = colorstr_to_idx(col_seg_str);
835
836 dinucs[i].second.fwd_string = idx;
837 dinucs[i].second.rev_string = rc_color_str(idx);
838 dinucs[i].second.first_in_string = s[1];
839 dinucs[i].second.last_in_string = s[half_splice_mer_len];
840 }
841 else
842 {
843 Dna5String seg_str = seqan::infixWithLength(ref_str,
844 pos + 2,
845 half_splice_mer_len);
846
847 stringstream ss(stringstream::in | stringstream::out);
848 string s;
849 ss << seg_str;
850 ss >> s;
851 uint64_t idx = dna5str_to_idx(s);
852
853 dinucs[i].second.fwd_string = idx;
854 dinucs[i].second.rev_string = rc_dna_str(idx);
855 }
856 }
857 }
858 };
859
860 struct PackedSplice
861 {
862 PackedSplice() : left(0u), seed(0u), right(0u) {}
863 PackedSplice(uint32_t l, uint64_t s, uint32_t r) : left(l), seed(s), right(r) {}
864 uint32_t left;
865 uint64_t seed;
866 uint32_t right;
867 };
868
869 // The second element of these pairs is the left (or right) side of a splice
870 // seed from a possible junction. The first element is the sequence flanking
871 // that seed
872 typedef pair<uint32_t, uint32_t> PackedSpliceHalf;
873
874 static inline std::string u32ToDna(uint32_t a, int len) {
875 char buf[17];
876 assert(len <= 16);
877 for(int i = 0; i < len; i++) {
878 buf[len-i-1] = "ACGT"[a & 3];
879 a >>= 2;
880 }
881 buf[len] = '\0';
882 return std::string(buf);
883 }
884
885
886 PackedSpliceHalf pack_left_splice_half(const string& seq,
887 uint32_t pos_in_l,
888 unsigned int seq_key_len)
889 {
890 const char* l = seq.c_str();
891 l += pos_in_l;
892
893 const char* left_end = l;
894 l -= 16;
895
896 assert (l + seq_key_len < seq.c_str() + seq.length());
897
898 PackedSpliceHalf packed_half = make_pair(0u,0u);
899
900 // Pack up to 32 bits (16 bases) of sequence into left
901 if (l < seq.c_str())
902 l = seq.c_str();
903 while (l < left_end)
904 {
905 packed_half.first <<= 2;
906 packed_half.first |= (0x3 & charToDna5[(size_t)*l]);
907 ++l;
908 }
909
910 // Pack up the seed bits
911 for (unsigned int i = 0; i < seq_key_len; ++i)
912 {
913 packed_half.second <<= 2;
914 packed_half.second |= (0x3 & charToDna5[(size_t)*(l + i)]);
915 }
916 return packed_half;
917 }
918
919
920 PackedSpliceHalf pack_right_splice_half(const string& seq,
921 uint32_t pos,
922 unsigned int seq_key_len)
923 {
924 const char* r = seq.c_str();
925 r += pos - seq_key_len;
926
927 PackedSpliceHalf packed_half = make_pair(0u,0u);
928
929 // Pack the seed bits
930 for (unsigned int i = 0; i < seq_key_len; ++i)
931 {
932 packed_half.second <<= 2;
933 packed_half.second |= (0x3 & charToDna5[(size_t)*(r + i)]);
934 }
935
936 r += seq_key_len;
937 // Now pack 32 bits (16 bases) of sequence into left
938 const char* right_end = r + 16;
939
940 if ((size_t)(right_end - seq.c_str()) > seq.length())
941 right_end = seq.c_str() + seq.length();
942 while (r < right_end)
943 {
944 packed_half.first <<= 2;
945 packed_half.first |= (0x3 & charToDna5[(size_t)*r]);
946 ++r;
947 }
948
949 return packed_half;
950 }
951
952 PackedSplice combine_splice_halves(const PackedSpliceHalf& left_half,
953 const PackedSpliceHalf& right_half,
954 int seq_key_len)
955 {
956 uint64_t seed = left_half.second << (seq_key_len << 1) | right_half.second;
957 return PackedSplice(left_half.first,seed, right_half.first);
958 }
959
960 PackedSplice pack_splice(const string& seq,
961 int l_pos_in_seq,
962 int r_pos_in_seq,
963 unsigned int seq_key_len)
964 {
965 const char* l = seq.c_str(); // l points to beginning of left exon sequence
966 l += l_pos_in_seq;
967
968 assert (l + seq_key_len < seq.c_str() + seq.length());
969
970 const char* r = seq.c_str(); // r points to beginning of right exon sequence
971 r += r_pos_in_seq - seq_key_len;
972 //r += 2; // r points to right side of junction;
973
974 uint64_t seed = 0;
975 uint64_t left = 0;
976 uint64_t right = 0;
977
978 // Pack up the seed bits
979 for (unsigned int i = 0; i < seq_key_len; ++i)
980 {
981 seed <<= 2;
982 seed |= (0x3 & charToDna5[(size_t)*(l + i)]);
983 }
984
985 for (unsigned int i = 0; i < seq_key_len; ++i)
986 {
987 seed <<= 2;
988 seed |= (0x3 & charToDna5[(size_t)*(r + i)]);
989 }
990
991 // Now pack 32 bits (16 bases) of sequence into left
992 const char* left_end = l;
993 l -= 16;
994 if (l < seq.c_str())
995 l = seq.c_str();
996 while (l < left_end)
997 {
998 left <<= 2;
999 left |= (0x3 & charToDna5[(size_t)*l]);
1000 ++l;
1001 }
1002
1003 r += seq_key_len;
1004 // Now pack 32 bits (16 bases) of sequence into left
1005 const char* right_end = r + 16;
1006
1007 if ((size_t)(right_end - seq.c_str()) > seq.length())
1008 right_end = seq.c_str() + seq.length();
1009 while (r < right_end)
1010 {
1011 right <<= 2;
1012 right |= (0x3 & charToDna5[(size_t)*r]);
1013 ++r;
1014 }
1015
1016 return PackedSplice((uint32_t)left,seed,(uint32_t)right);
1017 }
1018
1019
1020
1021 /* Represents a hit between a splice seed and a read. */
1022 // TODO: consider packing pos and meta into a single 32-bit int.
1023 struct ReadHit
1024 {
1025 ReadHit(uint32_t l, uint32_t r, uint32_t p, uint32_t m, bool rc)
1026 : left(l), right(r), pos(p), meta(m), reverse_complement(rc) {}
1027 uint32_t left; // 2-bits per base rep of the left remainder
1028 uint32_t right; //2-bits per base rep of the right remainder
1029 uint32_t pos; // position of the seed within the read
1030 uint32_t meta : 31;
1031 bool reverse_complement : 1;
1032 } __attribute__((packed));
1033
1034 // A MerTable maps k-mers to hits in indexed reads. See the comment for
1035 // mer_table
1036 typedef vector<ReadHit> ReadHitList;
1037 typedef vector<ReadHitList> MerTable;
1038
1039 size_t index_read(MerTable* mer_table,
1040 int seq_key_len,
1041 const string& seq,
1042 unsigned int read_num,
1043 bool reverse_complement,
1044 vector<uint64_t>& seeds)
1045 {
1046 // h is will hold the 2-bit-per-base representation of the k-mer seeds for
1047 // this read.
1048
1049 uint64_t seed = 0;
1050 bitset<256> left = 0;
1051 bitset<256> right = 0;
1052 const char* p = seq.c_str();
1053
1054 unsigned int seq_len = (int)seq.length();
1055 const char* seq_end = p + seq_len;
1056
1057 // Build the first seed
1058 while (p < seq.c_str() + (2 * seq_key_len))
1059 {
1060 seed <<= 2;
1061 seed |= (0x3 & charToDna5[(size_t)*p]);
1062 ++p;
1063 }
1064
1065 seeds.push_back(seed);
1066
1067 // Build the rest of them with a sliding window, adding successive bases
1068 // to the "right-remainder" word.
1069 while (p < seq_end)
1070 {
1071
1072 right <<= 2;
1073 right |= (0x3 & charToDna5[(size_t)*p]);
1074 ++p;
1075 }
1076
1077 size_t cap_increase = 0;
1078
1079 // At this point, seed contains the 5'-most 2*min_anchor_len bases of the
1080 // read, and right contains everthing else on the 3' end.
1081
1082 // This loop will construct successive seed, along with 32-bit words
1083 // containing the left and right remainders for each seed
1084 uint32_t i = 0;
1085 size_t new_hits = 0;
1086 do
1087 {
1088 // Let's not make an out-of-bounds write, if this fails the global
1089 // mer_table is too small
1090 assert (!mer_table || seed < mer_table->size());
1091
1092 // How many base pairs exist in the right remainder beyond what we have
1093 // space for ?
1094 int extra_right_bp = ((int)seq.length() - (i + 2 * seq_key_len)) - 16;
1095
1096 uint32_t hit_right;
1097 if (extra_right_bp > 0)
1098 {
1099 //bitset<32> tmp_hit_right = (right >> (extra_right_bp << 1));
1100 hit_right = (uint32_t)(right >> (extra_right_bp << 1)).to_ulong();
1101 }
1102 else
1103 {
1104 hit_right = (uint32_t)right.to_ulong();
1105 }
1106
1107 uint32_t hit_left = (uint32_t)((left << (256 - 32)) >> (256 - 32)).to_ulong();
1108
1109 if (mer_table)
1110 {
1111 size_t prev_cap = (*mer_table)[seed].capacity();
1112 (*mer_table)[seed].push_back(ReadHit(hit_left, hit_right,i, read_num, reverse_complement));
1113 cap_increase += ((*mer_table)[seed].capacity() - prev_cap) * sizeof (ReadHit);
1114 }
1115 new_hits++;
1116
1117 // Take the leftmost base of the seed and stick it into bp
1118 uint64_t bp = seed & (0x3uLL << ((seq_key_len << 2) - 2));
1119
1120 // Move that base down to the least significant bits of bp
1121 bp >>= ((seq_key_len << 2) - 2);
1122
1123 // And tack it onto the left remainder of the read
1124 left <<= 2;
1125 left |= bp;
1126
1127 // Now take the leftmost base of the right remainder and stick it into
1128 // the rightmost position of the seed
1129 uint32_t right_len = seq_len - (i + seq_key_len * 2);
1130 //bp = right & (0x3uLL << ((right_len - 1) << 1));
1131
1132 seed <<= 2;
1133 //cout << right << endl;
1134 bitset<256> tmp_right = (right >> ((right_len - 1) << 1));
1135 //cout <<tmp_right << endl;
1136 seed |= tmp_right.to_ulong();
1137 seed &= ~(0xFFFFFFFFFFFFFFFFuLL << (seq_key_len << 2));
1138
1139 seeds.push_back(seed);
1140
1141 //Now remove that leftmost base of the right remainder
1142 //right &= ~(0xFFFFFFFFFFFFFFFFuLL << (right_len - 1 << 1));
1143 if (right_len)
1144 {
1145 right.set(((right_len - 1) << 1), 0);
1146 right.set(((right_len - 1) << 1) + 1, 0);
1147 }
1148 ++i;
1149
1150 }while(i <= (size_t)(seq_end - seq.c_str()) - (2 * seq_key_len));
1151 return cap_increase;
1152 }
1153
1154
1155
1156 struct SeedExtension
1157 {
1158 SeedExtension(int lp,
1159 int rp,
1160 int rd_p,
1161 int le,
1162 int re,
1163 int mm) :
1164 l_pos_in_ref(lp),
1165 r_pos_in_ref(rp),
1166 read_pos(rd_p),
1167 left_extent(le),
1168 right_extent(re),
1169 mismatches(mm){}
1170
1171 int l_pos_in_ref;
1172 int r_pos_in_ref;
1173 int read_pos;
1174 int left_extent;
1175 int right_extent;
1176 int mismatches;
1177 };
1178
1179 pair<string::size_type, int> left_extend(const string& ref,
1180 const string& read,
1181 int ref_pos,
1182 int read_pos,
1183 int num_mismatches)
1184 {
1185 string::size_type ext = 0;
1186 int mm_encountered = 0;
1187 while(ref_pos >= 0 &&
1188 read_pos >= 0)
1189 {
1190 //char ref_char = ref[ref_pos];
1191 //char read_char = read[read_pos];
1192 if (ref[ref_pos] != read[read_pos])
1193 {
1194 if (mm_encountered + 1 > num_mismatches)
1195 return make_pair(ext, mm_encountered);
1196 mm_encountered++;
1197 }
1198 ext++;
1199
1200 --ref_pos;
1201 --read_pos;
1202 }
1203
1204 return make_pair(ext, mm_encountered);
1205 }
1206
1207 pair<string::size_type, int> right_extend(const string& ref,
1208 const string& read,
1209 int ref_pos,
1210 int read_pos,
1211 int num_mismatches)
1212 {
1213 string::size_type ext = 0;
1214 int mm_encountered = 0;
1215 while(ref_pos < (int)ref.size() &&
1216 read_pos < (int)read.size())
1217 {
1218 if (ref[ref_pos] != read[read_pos])
1219 {
1220 if (mm_encountered + 1 > num_mismatches)
1221 return make_pair(ext, mm_encountered);
1222 mm_encountered++;
1223 }
1224
1225 ext++;
1226
1227 ++ref_pos;
1228 ++read_pos;
1229 }
1230
1231 return make_pair(ext, mm_encountered);
1232 }
1233
1234
1235 void extend_from_seeds(vector<SeedExtension>& extensions,
1236 const PackedSplice& p,
1237 const MerTable& mer_table,
1238 const string& ref,
1239 const string& read,
1240 size_t l_pos_in_ref,
1241 size_t r_pos_in_ref,
1242 int seq_key_len)
1243 {
1244 assert(p.seed < mer_table.size());
1245 const ReadHitList& hl = mer_table[p.seed];
1246
1247 for (size_t hit = 0; hit < hl.size(); ++hit)
1248 {
1249 const ReadHit& rh = hl[hit];
1250 uint32_t pos = rh.pos;
1251
1252 pair<string::size_type, int> left_extension;
1253 pair<string::size_type, int> right_extension;
1254 left_extension = left_extend(ref, read, l_pos_in_ref - seq_key_len + 1, pos, 2);
1255 right_extension = right_extend(ref, read, r_pos_in_ref + seq_key_len, pos + 2 * seq_key_len, 2);
1256 extensions.push_back(SeedExtension(l_pos_in_ref,
1257 r_pos_in_ref,
1258 pos + seq_key_len,
1259 left_extension.first,
1260 right_extension.first,
1261 left_extension.second + right_extension.second));
1262
1263 }
1264 }
1265
1266 typedef pair<size_t, PackedSpliceHalf> SpliceHalf;
1267
1268 void get_seed_extensions(const string& ref,
1269 const string& read,
1270 int seq_key_len,
1271 MerTable& mer_table,
1272 vector<SpliceHalf>& donors,
1273 vector<SpliceHalf>& acceptors,
1274 vector<SeedExtension>& extensions)
1275 {
1276 for (size_t d = 0; d < donors.size(); ++d)
1277 {
1278 bool broke_out = false;
1279
1280 // start pos is a lower bound on downstream acceptor positions
1281 // to consider
1282 size_t start_pos = donors[d].first + min_report_intron_length;
1283
1284 SpliceHalf dummy = make_pair(start_pos,PackedSpliceHalf());
1285 vector<SpliceHalf>::iterator lb = upper_bound(acceptors.begin(),
1286 acceptors.end(),
1287 dummy);
1288
1289 if (lb == acceptors.end())
1290 break;
1291 for (size_t a = lb - acceptors.begin();
1292 a < acceptors.size();
1293 ++a)
1294 {
1295 if (acceptors[a].first - donors[d].first > (size_t)max_microexon_stretch)
1296 {
1297 broke_out = true;
1298 break;
1299 }
1300
1301 size_t l_pos_in_ref = donors[d].first - 1;
1302 size_t r_pos_in_ref = acceptors[a].first + 2;
1303 PackedSplice p = combine_splice_halves(donors[d].second,
1304 acceptors[a].second,
1305 seq_key_len);
1306 extend_from_seeds(extensions,
1307 p,
1308 mer_table,
1309 ref,
1310 read,
1311 l_pos_in_ref,
1312 r_pos_in_ref,
1313 seq_key_len);
1314 }
1315 if (broke_out)
1316 continue;
1317 }
1318
1319 }
1320
1321 void hits_from_seed_extension(uint32_t ref_id,
1322 int ref_offset,
1323 uint64_t insert_id,
1324 bool antisense,
1325 vector<SeedExtension>& extensions,
1326 vector<BowtieHit>& hits_out,
1327 int left_read_edge,
1328 int right_read_edge,
1329 int seq_key_len)
1330 {
1331 for (size_t i = 0; i < extensions.size(); ++i)
1332 {
1333 SeedExtension& s = extensions[i];
1334 if (s.read_pos >= right_read_edge ||
1335 s.read_pos < left_read_edge)
1336 continue;
1337 if (s.read_pos - seq_key_len - s.left_extent <= left_read_edge &&
1338 s.read_pos + seq_key_len + s.right_extent >= right_read_edge && s.mismatches <= 2 )
1339 {
1340 vector<CigarOp> cigar;
1341 int off_adjust;
1342 if (antisense)
1343 {
1344 CigarOp m1 = CigarOp(MATCH, s.read_pos - left_read_edge);
1345 CigarOp skip = CigarOp(REF_SKIP, s.r_pos_in_ref - s.l_pos_in_ref);
1346 CigarOp m2 = CigarOp(MATCH, right_read_edge - s.read_pos);
1347 cigar.push_back(m1);
1348 cigar.push_back(skip);
1349 cigar.push_back(m2);
1350 off_adjust = m1.length;
1351 }
1352 else
1353 {
1354 CigarOp m1 = CigarOp(MATCH, s.read_pos - left_read_edge + 1);
1355 CigarOp skip = CigarOp(REF_SKIP, s.r_pos_in_ref - s.l_pos_in_ref);
1356 CigarOp m2 = CigarOp(MATCH, right_read_edge - s.read_pos - 1);
1357 cigar.push_back(m1);
1358 cigar.push_back(skip);
1359 cigar.push_back(m2);
1360 off_adjust = m1.length;
1361 }
1362 // daehwan - check this
1363 bool end = false;
1364 BowtieHit bh(ref_id,
1365 insert_id,
1366 ref_offset + s.l_pos_in_ref - off_adjust + 1,
1367 cigar,
1368 antisense,
1369 false,
1370 s.mismatches,
1371 0,
1372 end);
1373 hits_out.push_back(bh);
1374 }
1375 }
1376 }
1377
1378 void align(uint32_t ref_id,
1379 uint64_t insert_id,
1380 bool antisense,
1381 const string& ref,
1382 const string& read,
1383 int ref_offset,
1384 int left_read_edge,
1385 int right_read_edge,
1386 MerTable& mer_table,
1387 int seq_key_len,
1388 vector<BowtieHit>& hits_out)
1389 {
1390 // Reserve an entry for each k-mer we might see
1391 size_t mer_table_size = 1 << ((seq_key_len << 1)<<1);
1392
1393 mer_table.resize(mer_table_size);
1394
1395 vector<uint64_t> seeds;
1396 index_read(&mer_table, seq_key_len, read, 0, antisense, seeds);
1397
1398 vector<SpliceHalf> forward_donors;
1399 vector<SpliceHalf> forward_acceptors;
1400 vector<SpliceHalf> reverse_donors;
1401 vector<SpliceHalf> reverse_acceptors;
1402
1403 const string& seq = ref;
1404
1405 unsigned int pos = 0;
1406
1407
1408 for (size_t z = seq_key_len + 1; z < seq.length() - seq_key_len - 2; ++z)
1409 {
1410 char l = seq[z - 1];
1411 char r = seq[z];
1412 if (l == 'G' && r == 'T')
1413 {
1414 size_t donor_pos = pos + z - 1;
1415 size_t s = donor_pos - seq_key_len;
1416 PackedSpliceHalf p = pack_left_splice_half(seq, s, seq_key_len);
1417 forward_donors.push_back(make_pair(donor_pos,p));
1418 }
1419 if (l == 'A' && r == 'G')
1420 {
1421 size_t acceptor_pos = pos + z - 1;
1422 size_t s = acceptor_pos + 2 + seq_key_len;
1423 PackedSpliceHalf p = pack_right_splice_half(seq, s, seq_key_len);
1424 forward_acceptors.push_back(make_pair(acceptor_pos,p));
1425 }
1426 if (l == 'C' && r == 'T')
1427 {
1428 size_t acceptor_pos = pos + z - 1;
1429 size_t s = acceptor_pos - seq_key_len;
1430 PackedSpliceHalf p = pack_left_splice_half(seq, s, seq_key_len);
1431 reverse_acceptors.push_back(make_pair(pos + z - 1,p));
1432 }
1433 if (l == 'A' && r == 'C')
1434 {
1435 size_t donor_pos = pos + z - 1;
1436 size_t s = donor_pos + 2 + seq_key_len;
1437 PackedSpliceHalf p = pack_right_splice_half(seq, s, seq_key_len);
1438 reverse_donors.push_back(make_pair(donor_pos,p));
1439 }
1440 }
1441
1442 vector<SeedExtension> fwd_extensions;
1443 get_seed_extensions(seq,
1444 read,
1445 seq_key_len,
1446 mer_table,
1447 forward_donors,
1448 forward_acceptors,
1449 fwd_extensions);
1450
1451 hits_from_seed_extension(ref_id,
1452 ref_offset,
1453 insert_id,
1454 antisense,
1455 fwd_extensions,
1456 hits_out,
1457 left_read_edge,
1458 right_read_edge,
1459 seq_key_len);
1460
1461 //fprintf(stderr, "Found %d seed hits\n", fwd_extensions.size());
1462
1463 vector<SeedExtension> rev_extensions;
1464 get_seed_extensions(seq,
1465 read,
1466 seq_key_len,
1467 mer_table,
1468 reverse_donors,
1469 reverse_acceptors,
1470 rev_extensions);
1471
1472 hits_from_seed_extension(ref_id,
1473 ref_offset,
1474 insert_id,
1475 antisense,
1476 rev_extensions,
1477 hits_out,
1478 left_read_edge,
1479 right_read_edge,
1480 seq_key_len);
1481
1482 for (size_t i = 0; i < seeds.size(); ++i)
1483 mer_table[seeds[i]].clear();
1484 }
1485
1486 int extension_mismatches = 0;
1487
1488
1489 bool left_extendable_junction(uint64_t upstream_dna_str,
1490 size_t key,
1491 size_t splice_mer_len,
1492 size_t min_ext_len)
1493 {
1494 vector<MerExtension>& exts = extensions[key];
1495 for (size_t i = 0; i < exts.size(); ++i)
1496 {
1497 const MerExtension& ext = exts[i];
1498 if (ext.left_ext_len < min_ext_len)
1499 continue;
1500 uint64_t upstream = upstream_dna_str & ~(0xFFFFFFFFFFFFFFFFull << (ext.left_ext_len << 1));
1501 int mism = mismatching_bases(ext.left_dna_str, upstream, ext.left_ext_len, extension_mismatches);
1502 if (mism <= extension_mismatches)
1503 return true;
1504 }
1505 return false;
1506 }
1507
1508 bool right_extendable_junction(uint64_t downstream_dna_str,
1509 size_t key,
1510 size_t splice_mer_len,
1511 size_t min_ext_len)
1512 {
1513 vector<MerExtension>& exts = extensions[key];
1514 for (size_t i = 0; i < exts.size(); ++i)
1515 {
1516 const MerExtension& ext = exts[i];
1517 if (ext.right_ext_len < min_ext_len)
1518 continue;
1519 uint64_t mask = ~(0xFFFFFFFFFFFFFFFFull >> (ext.right_ext_len << 1));
1520 uint64_t downstream = downstream_dna_str & mask;
1521 downstream >>= ((32 - ext.right_ext_len) << 1);
1522 int mism = mismatching_bases(ext.right_dna_str, downstream, ext.right_ext_len, extension_mismatches);
1523
1524 if (mism <= extension_mismatches)
1525 return true;
1526 }
1527 return false;
1528 }
1529
1530 uint32_t junction_key(uint64_t upstream_dna_str,
1531 uint64_t downstream_dna_str,
1532 size_t splice_mer_len)
1533 {
1534 uint64_t upstream_mask = ~(0xFFFFFFFFFFFFFFFFull << (splice_mer_len));
1535 uint64_t upstream_key_half = upstream_dna_str & upstream_mask;
1536 uint64_t downstream_mask = ~(0xFFFFFFFFFFFFFFFFull >> (splice_mer_len));
1537 uint64_t downstream_key_half = (downstream_dna_str & downstream_mask) >> (64 - splice_mer_len);
1538 uint32_t key = ((uint32_t)upstream_key_half << splice_mer_len) | (uint32_t)downstream_key_half;
1539 return key;
1540 }
1541
1542 bool extendable_junction(uint64_t upstream_dna_str,
1543 uint64_t downstream_dna_str,
1544 size_t splice_mer_len,
1545 size_t min_ext_len,
1546 bool reverse,
1547 char last_in_upstream = 'N',
1548 char first_in_downstream = 'N')
1549 {
1550 if (color)
1551 {
1552 string two_bp; two_bp.push_back(last_in_upstream); two_bp.push_back(first_in_downstream);
1553 string color = convert_bp_to_color(two_bp, true);
1554 char num = (color[0] - '0') & 0x3;
1555
1556 if (reverse)
1557 {
1558 upstream_dna_str = (upstream_dna_str >> 2) << 2;
1559 upstream_dna_str |= (uint64_t)num;
1560 }
1561 else
1562 {
1563 downstream_dna_str = (downstream_dna_str << 2) >> 2;
1564 downstream_dna_str |= ((uint64_t)num << 62);
1565 }
1566 }
1567
1568 uint32_t key = junction_key(upstream_dna_str,
1569 downstream_dna_str,
1570 splice_mer_len);
1571
1572 upstream_dna_str >>= splice_mer_len;
1573 downstream_dna_str <<= splice_mer_len;
1574
1575 bool extendable = (left_extendable_junction(upstream_dna_str,
1576 key, splice_mer_len, min_ext_len) ||
1577 right_extendable_junction(downstream_dna_str,
1578 key, splice_mer_len, min_ext_len));
1579 return extendable;
1580 }
1581
1582 typedef std::set<Junction, skip_count_lt> PotentialJuncs;
1583
1584 struct RecordExtendableJuncs
1585 {
1586 void record(uint32_t ref_id,
1587 const vector<pair<size_t, DnaSpliceStrings> >& left_sites,
1588 const vector<pair<size_t, DnaSpliceStrings> >& right_sites,
1589 bool antisense,
1590 PotentialJuncs& juncs,
1591 int min_intron,
1592 int max_intron,
1593 size_t max_juncs,
1594 size_t half_splice_mer_len)
1595 {
1596
1597 size_t splice_mer_len = 2 * half_splice_mer_len;
1598
1599 size_t curr_R = 0;
1600 for (size_t L = 0; L < left_sites.size(); ++L)
1601 {
1602 while (curr_R < right_sites.size() &&
1603 right_sites[curr_R].first < left_sites[L].first + min_intron)
1604 {
1605 curr_R++;
1606 }
1607
1608 size_t left_pos = left_sites[L].first;
1609 size_t max_right_pos = left_pos + max_intron;
1610 for (size_t R = curr_R;
1611 R < right_sites.size() && right_sites[R].first < max_right_pos; ++R)
1612 {
1613 uint64_t upstream_dna_str = left_sites[L].second.fwd_string;
1614 char last_in_upstream = left_sites[L].second.last_in_string;
1615 uint64_t downstream_dna_str = right_sites[R].second.fwd_string;
1616 char first_in_downstream = right_sites[R].second.first_in_string;
1617 uint64_t rc_upstream_dna_str = left_sites[L].second.rev_string;
1618 uint64_t rc_downstream_dna_str = right_sites[R].second.rev_string;
1619
1620 if (extendable_junction(upstream_dna_str,
1621 downstream_dna_str, splice_mer_len, 7, false,
1622 last_in_upstream, first_in_downstream) ||
1623 extendable_junction(rc_downstream_dna_str,
1624 rc_upstream_dna_str, splice_mer_len, 7, true,
1625 last_in_upstream, first_in_downstream))
1626 {
1627 juncs.insert(Junction(ref_id,
1628 left_sites[L].first - 1,
1629 right_sites[R].first + 2,
1630 antisense,
1631 R - curr_R));
1632 }
1633 if (juncs.size() > max_juncs)
1634 juncs.erase(*(juncs.rbegin()));
1635 }
1636 }
1637 }
1638 };
1639
1640 struct RecordAllJuncs
1641 {
1642 void record(uint32_t ref_id,
1643 const vector<pair<size_t, DnaSpliceStrings> >& left_sites,
1644 const vector<pair<size_t, DnaSpliceStrings> >& right_sites,
1645 bool antisense,
1646 PotentialJuncs& juncs,
1647 int min_intron,
1648 int max_intron,
1649 size_t max_juncs,
1650 size_t half_splice_mer_len)
1651 {
1652 size_t curr_R = 0;
1653 for (size_t L = 0; L < left_sites.size(); ++L)
1654 {
1655 while (curr_R < right_sites.size() &&
1656 right_sites[curr_R].first < left_sites[L].first + min_intron)
1657 {
1658 curr_R++;
1659 }
1660
1661 size_t left_pos = left_sites[L].first;
1662 size_t max_right_pos = left_pos + max_intron;
1663 for (size_t R = curr_R;
1664 R < right_sites.size() && right_sites[R].first < max_right_pos; ++R)
1665 {
1666 Junction j(ref_id,
1667 left_sites[L].first - 1,
1668 right_sites[R].first + 2,
1669 antisense,
1670 R - curr_R);
1671
1672 juncs.insert(j);
1673
1674 if (juncs.size() > max_juncs)
1675 juncs.erase(*(juncs.rbegin()));
1676 }
1677 }
1678 }
1679 };
1680
1681 struct RecordSegmentJuncs
1682 {
1683 void record(uint32_t ref_id,
1684 const vector<pair<size_t, DnaSpliceStrings> >& left_sites,
1685 const vector<pair<size_t, DnaSpliceStrings> >& right_sites,
1686 bool antisense,
1687 PotentialJuncs& juncs,
1688 int min_intron,
1689 int max_intron,
1690 size_t max_juncs,
1691 size_t half_splice_mer_len)
1692 {
1693 if (left_sites.size() != right_sites.size())
1694 return;
1695
1696 for (size_t i = 0; i < left_sites.size(); ++i)
1697 {
1698 Junction j(ref_id,
1699 left_sites[i].first - 1,
1700 right_sites[i].first + 2,
1701 antisense);
1702
1703 juncs.insert(j);
1704
1705 if (juncs.size() > max_juncs)
1706 juncs.erase(*(juncs.rbegin()));
1707 }
1708 }
1709 };
1710
1711 struct ButterflyKey
1712 {
1713 uint32_t pos;
1714 uint32_t key;
1715
1716 ButterflyKey(uint32_t p, uint32_t k) : pos(p), key(k) {}
1717
1718 bool operator<(const ButterflyKey& rhs) const
1719 {
1720 if (key != rhs.key)
1721 return key < rhs.key;
1722 if (pos != rhs.pos)
1723 return pos < rhs.pos;
1724 return false;
1725 }
1726
1727 bool operator==(const ButterflyKey& rhs) const
1728 {
1729 return pos == rhs.pos && key == rhs.key;
1730 }
1731 };
1732
1733 uint32_t get_left_butterfly_key(uint64_t upstream_key,
1734 const MerExtension& ext,
1735 size_t half_splice_mer_len)
1736 {
1737 uint64_t key = ext.right_dna_str >> ((ext.right_ext_len - half_splice_mer_len) << 1);
1738 uint64_t mask = ~(0xFFFFFFFFFFFFFFFFull << (half_splice_mer_len << 1));
1739 uint64_t top_half = upstream_key & mask;
1740 key |= (top_half << (half_splice_mer_len << 1));
1741 return (uint32_t)key;
1742 }
1743
1744 uint32_t get_right_butterfly_key(uint64_t downstream_key,
1745 const MerExtension& ext,
1746 size_t half_splice_mer_len)
1747 {
1748 uint64_t mask = ~(0xFFFFFFFFFFFFFFFFull << (half_splice_mer_len << 1));
1749 uint64_t key = (ext.left_dna_str & mask) << (half_splice_mer_len << 1);
1750 uint64_t bottom_half = (downstream_key >> (half_splice_mer_len << 1));
1751 key |= bottom_half;
1752 return (uint32_t)key;
1753 }
1754
1755 struct RecordButterflyJuncs
1756 {
1757 void record(uint32_t ref_id,
1758 const vector<pair<size_t, DnaSpliceStrings> >& all_left_sites,
1759 const vector<pair<size_t, DnaSpliceStrings> >& all_right_sites,
1760 bool antisense,
1761 PotentialJuncs& juncs,
1762 int min_intron,
1763 int max_intron,
1764 size_t max_juncs,
1765 size_t half_splice_mer_len)
1766 {
1767 size_t key_length = 2 * half_splice_mer_len;
1768 size_t extension_length = butterfly_overhang;
1769 uint64_t bottom_bit_mask = ~(0xFFFFFFFFFFFFFFFFull << (key_length<<1));
1770 uint64_t top_bit_mask = ~(0xFFFFFFFFFFFFFFFFull >> (key_length<<1));
1771
1772 if (all_left_sites.empty() || all_right_sites.empty())
1773 return;
1774
1775 size_t last_site = max(all_left_sites.back().first,
1776 all_right_sites.back().first);
1777
1778 size_t curr_left_site = 0;
1779 size_t curr_right_site = 0;
1780 for (size_t window_left_edge = 0;
1781 window_left_edge < last_site;
1782 window_left_edge += max_intron)
1783 {
1784 //fprintf(stderr, "\twindow %lu - %lu\n", window_left_edge, window_left_edge + 2 * max_intron);
1785 vector<pair<size_t, DnaSpliceStrings> > left_sites;
1786 vector<pair<size_t, DnaSpliceStrings> > right_sites;
1787
1788 while(curr_left_site < all_left_sites.size() &&
1789 all_left_sites[curr_left_site].first < window_left_edge)
1790 {
1791 curr_left_site++;
1792 }
1793
1794 while(curr_right_site < all_right_sites.size() &&
1795 all_right_sites[curr_right_site].first < window_left_edge)
1796 {
1797 curr_right_site++;
1798 }
1799
1800 for (size_t ls = curr_left_site; ls < all_left_sites.size(); ++ls)
1801 {
1802 if (all_left_sites[ls].first < window_left_edge + 2 * max_intron)
1803 {
1804 left_sites.push_back(all_left_sites[ls]);
1805 }
1806 }
1807
1808 for (size_t rs = curr_right_site; rs < all_right_sites.size(); ++rs)
1809 {
1810 if (all_right_sites[rs].first < window_left_edge + 2 * max_intron)
1811 {
1812 right_sites.push_back(all_right_sites[rs]);
1813 }
1814 }
1815
1816 vector<ButterflyKey> left_keys;
1817 for (size_t L = 0; L < left_sites.size(); ++L)
1818 {
1819 uint64_t fwd_upstream_dna_str = left_sites[L].second.fwd_string;
1820 uint64_t fwd_upstream_key = fwd_upstream_dna_str & bottom_bit_mask;
1821
1822 assert (fwd_upstream_key < extensions.size());
1823
1824 vector<MerExtension>& fwd_exts = extensions[fwd_upstream_key];
1825 for (size_t i = 0; i < fwd_exts.size(); ++i)
1826 {
1827 const MerExtension& ext = fwd_exts[i];
1828 if (ext.right_ext_len < extension_length)
1829 continue;
1830
1831 /*
1832 < f_u_key ><ext.right>
1833 NNNNNNNNNN GT
1834 */
1835
1836 // take the top bits of the right extension
1837 uint64_t key = ext.right_dna_str >> ((ext.right_ext_len - extension_length) << 1);
1838
1839 // and the bottom bits of the site key
1840 uint64_t mask = ~(0xFFFFFFFFFFFFFFFFull << (extension_length << 1));
1841 uint64_t top_half = fwd_upstream_key & mask;
1842
1843 // and cat them together
1844 key |= (top_half << (extension_length << 1));
1845 left_keys.push_back(ButterflyKey((uint32_t)left_sites[L].first, key));
1846 }
1847
1848 uint64_t rev_upstream_dna_str = left_sites[L].second.rev_string;
1849 uint64_t rev_upstream_key = (rev_upstream_dna_str & top_bit_mask) >> (64 - (key_length<<1));
1850
1851 assert (rev_upstream_key < extensions.size());
1852
1853 vector<MerExtension>& rev_exts = extensions[rev_upstream_key];
1854 for (size_t i = 0; i < rev_exts.size(); ++i)
1855 {
1856 const MerExtension& ext = rev_exts[i];
1857 if (ext.left_ext_len < extension_length)
1858 continue;
1859
1860 /*
1861 < r_u_key ><ext.left>
1862 NNNNNNNNNN GT
1863 */
1864
1865 // reverse complement the left extension, and we will need
1866 // what were the bottom bits. these become the top bits in the
1867 // rc.
1868 uint64_t ext_str = color ? rc_color_str(ext.left_dna_str) : rc_dna_str(ext.left_dna_str);
1869 ext_str >>= 64 - (ext.left_ext_len << 1);
1870
1871 // now take the top bits of the rc, make them the bottom of
1872 // the key
1873 uint64_t key = ext_str >> ((ext.left_ext_len - extension_length) << 1);
1874
1875 // now add in the seed key bottom bits, making them the top of
1876 // the key
1877 uint64_t mask = ~(0xFFFFFFFFFFFFFFFFull << (extension_length << 1));
1878 uint64_t top_half = fwd_upstream_key & mask;
1879 key |= (top_half << (extension_length << 1));
1880 left_keys.push_back(ButterflyKey((uint32_t)left_sites[L].first, key));
1881 }
1882 }
1883 sort (left_keys.begin(), left_keys.end());
1884 vector<ButterflyKey>::iterator new_end = unique(left_keys.begin(), left_keys.end());
1885 left_keys.erase(new_end, left_keys.end());
1886
1887 vector<ButterflyKey> right_keys;
1888 for (size_t R = 0; R < right_sites.size(); ++R)
1889 {
1890 uint64_t fwd_downstream_dna_str = right_sites[R].second.fwd_string;
1891 uint64_t fwd_downstream_key = (fwd_downstream_dna_str & top_bit_mask) >> (64 - (key_length<<1));
1892
1893 assert (fwd_downstream_key < extensions.size());
1894
1895 vector<uint64_t> fwd_downstream_keys;
1896 if (color)
1897 {
1898 for(size_t color_value = 0; color_value < 4; ++color_value)
1899 {
1900 uint64_t tmp_key = (fwd_downstream_key << 2) >> 2 | (color_value << ((key_length - 1) << 1));
1901 fwd_downstream_keys.push_back(tmp_key);
1902 }
1903 }
1904 else
1905 {
1906 fwd_downstream_keys.push_back(fwd_downstream_key);
1907 }
1908
1909 for(size_t key = 0; key < fwd_downstream_keys.size(); ++key)
1910 {
1911 uint64_t tmp_fwd_downstream_key = fwd_downstream_keys[key];
1912 vector<MerExtension>& fwd_exts = extensions[tmp_fwd_downstream_key];
1913 for (size_t i = 0; i < fwd_exts.size(); ++i)
1914 {
1915 const MerExtension& ext = fwd_exts[i];
1916 if (ext.left_ext_len < extension_length)
1917 continue;
1918
1919 /*
1920 <ext.left>< f_d_key >
1921 AG NNNNNNNNNN
1922 */
1923
1924 // take the bottom bits of the left extension, making them the
1925 // top of the key.
1926 uint64_t mask = ~(0xFFFFFFFFFFFFFFFFull << (extension_length << 1));
1927 uint64_t key = (ext.left_dna_str & mask) << (extension_length << 1);
1928
1929 // add in the top bits of the seed key, making them the bottom bits
1930 // of the key.
1931 uint64_t bottom_half = (tmp_fwd_downstream_key >> ((key_length - extension_length) << 1));
1932 key |= bottom_half;
1933 right_keys.push_back(ButterflyKey((uint32_t)right_sites[R].first, key));
1934 }
1935 }
1936
1937 uint64_t rev_downstream_dna_str = right_sites[R].second.rev_string;
1938 uint64_t rev_downstream_key = rev_downstream_dna_str & bottom_bit_mask;
1939
1940 assert (rev_downstream_key < extensions.size());
1941
1942 vector<uint64_t> rev_downstream_keys;
1943 if (color)
1944 {
1945 for(size_t color_value = 0; color_value < 4; ++color_value)
1946 {
1947 uint64_t tmp_key = (rev_downstream_key >> 2) << 2 | color_value;
1948 rev_downstream_keys.push_back(tmp_key);
1949 }
1950 }
1951 else
1952 {
1953 rev_downstream_keys.push_back(rev_downstream_key);
1954 }
1955
1956 for(size_t key = 0; key < rev_downstream_keys.size(); ++key)
1957 {
1958 uint64_t tmp_rev_downstream_key = rev_downstream_keys[key];
1959 uint64_t tmp_fwd_downstream_key = fwd_downstream_key;
1960 if (color)
1961 {
1962 tmp_fwd_downstream_key = rc_color_str(tmp_rev_downstream_key) >> (64 - (key_length << 1));
1963 }
1964
1965 vector<MerExtension>& rev_exts = extensions[tmp_rev_downstream_key];
1966 for (size_t i = 0; i < rev_exts.size(); ++i)
1967 {
1968 const MerExtension& ext = rev_exts[i];
1969 if (ext.right_ext_len < extension_length)
1970 continue;
1971
1972 /*
1973 <ext.right>< r_d_key >
1974 AG NNNNNNNNNN
1975 */
1976
1977 // reverse complement the right_extension. we want the
1978 // top bits of the extension, but these become the bottom bits
1979 // under the rc.
1980 uint64_t ext_str = color ? rc_color_str(ext.right_dna_str) : rc_dna_str(ext.right_dna_str);
1981 ext_str >>= 64 - (ext.right_ext_len << 1);
1982
1983 // take the bottom bits of the rc and make it the top of the key
1984 uint64_t key = ext_str << (extension_length << 1);
1985
1986 // take the top bits of the seed key and make them the bottom
1987 // of the key.
1988 uint64_t bottom_half = (tmp_fwd_downstream_key >> ((key_length - extension_length) << 1));
1989 key |= bottom_half;
1990
1991 right_keys.push_back(ButterflyKey((uint32_t)right_sites[R].first, key));
1992 }
1993 }
1994 }
1995
1996 sort (right_keys.begin(), right_keys.end());
1997 new_end = unique(right_keys.begin(), right_keys.end());
1998 right_keys.erase(new_end, right_keys.end());
1999
2000 size_t lk = 0;
2001 size_t rk = 0;
2002
2003 while (lk < left_keys.size() && rk < right_keys.size())
2004 {
2005 while (lk < left_keys.size() &&
2006 left_keys[lk].key < right_keys[rk].key) { ++lk; }
2007
2008 if (lk == left_keys.size())
2009 break;
2010
2011 while (rk < right_keys.size() &&
2012 right_keys[rk].key < left_keys[lk].key) { ++rk; }
2013
2014 if (rk == right_keys.size())
2015 break;
2016
2017 if (lk < left_keys.size() && rk < right_keys.size() &&
2018 right_keys[rk].key == left_keys[lk].key)
2019 {
2020
2021 size_t k = right_keys[rk].key;
2022 size_t lk_end = lk;
2023 size_t rk_end = rk;
2024 while (rk_end < right_keys.size() && right_keys[rk_end].key == k) {++rk_end;}
2025 while (lk_end < left_keys.size() && left_keys[lk_end].key == k) {++lk_end;}
2026
2027 size_t tmp_lk = lk;
2028
2029 while (tmp_lk < lk_end)
2030 {
2031 size_t tmp_rk = rk;
2032 while (tmp_rk < rk_end)
2033 {
2034 int donor = (int)left_keys[tmp_lk].pos - 1;
2035 int acceptor = (int)right_keys[tmp_rk].pos + 2;
2036
2037 if (acceptor - donor > min_intron && acceptor - donor < max_intron)
2038 {
2039 Junction j(ref_id,
2040 donor,
2041 acceptor,
2042 antisense,
2043 acceptor - donor); // just prefer shorter introns
2044 juncs.insert(j);
2045 if (juncs.size() > max_juncs)
2046 {
2047 juncs.erase(*(juncs.rbegin()));
2048 }
2049 }
2050 ++tmp_rk;
2051 }
2052
2053 ++tmp_lk;
2054 }
2055
2056 lk = lk_end;
2057 rk = rk_end;
2058 }
2059 }
2060 }
2061 }
2062 };
2063
2064 template <class JunctionRecorder>
2065 void juncs_from_ref_segs(RefSequenceTable& rt,
2066 vector<RefSeg>& expected_don_acc_windows,
2067 PotentialJuncs& juncs,
2068 const DnaString& donor_dinuc,
2069 const DnaString& acceptor_dinuc,
2070 int max_intron,
2071 int min_intron,
2072 size_t max_juncs,
2073 bool talkative,
2074 size_t half_splice_mer_len)
2075 {
2076
2077 typedef map<uint32_t, IntronMotifs> MotifMap;
2078
2079 MotifMap ims;
2080
2081 seqan::DnaStringReverseComplement rev_donor_dinuc(donor_dinuc);
2082 seqan::DnaStringReverseComplement rev_acceptor_dinuc(acceptor_dinuc);
2083
2084 if (talkative)
2085 fprintf(stderr, "Collecting potential splice sites in islands\n");
2086
2087 //
2088 bool all_both = true;
2089 for (size_t r = 0; r < expected_don_acc_windows.size(); ++r)
2090 {
2091 const RefSeg& seg = expected_don_acc_windows[r];
2092
2093 if (seg.points_where != POINT_DIR_BOTH)
2094 all_both = false;
2095
2096 RefSequenceTable::Sequence* ref_str = rt.get_seq(seg.ref_id);
2097
2098 if (!ref_str)
2099 continue;
2100
2101 bool skip_fwd = false;
2102 bool skip_rev = false;
2103
2104 if (library_type == FR_FIRSTSTRAND)
2105 {
2106 if (seg.read == READ_LEFT)
2107 {
2108 if (seg.antisense) skip_rev = true;
2109 else skip_fwd = true;
2110 }
2111 else if(seg.read == READ_RIGHT)
2112 {
2113 if (seg.antisense) skip_fwd = true;
2114 else skip_rev = true;
2115 }
2116 }
2117
2118 if (library_type == FR_SECONDSTRAND)
2119 {
2120 if (seg.read == READ_LEFT)
2121 {
2122 if (seg.antisense) skip_fwd = true;
2123 else skip_rev = true;
2124 }
2125 else if(seg.read == READ_RIGHT)
2126 {
2127 if (seg.antisense) skip_rev = true;
2128 else skip_fwd = true;
2129 }
2130 }
2131
2132 pair<map<uint32_t, IntronMotifs>::iterator, bool> ret =
2133 ims.insert(make_pair(seg.ref_id, IntronMotifs(seg.ref_id)));
2134
2135 IntronMotifs& motifs = ret.first->second;
2136
2137 int left_color_offset = 0, right_color_offset = 0;
2138 if (color)
2139 {
2140 if (seg.antisense)
2141 right_color_offset = 1;
2142 else
2143 left_color_offset = -1;
2144 }
2145
2146 if (seg.left + left_color_offset < 0 || seg.right + right_color_offset >= (int)length(*ref_str) - 1)
2147 continue;
2148
2149 DnaString org_seg_str = seqan::infix(*ref_str, seg.left + left_color_offset, seg.right + right_color_offset);
2150
2151 String<char> seg_str;
2152 assign(seg_str, org_seg_str);
2153
2154 // daehwan
2155 if (bDebug)
2156 {
2157 cout << "coord: " << seg.left << " " << seg.right << endl
2158 << "seg_str: " << seg_str << endl;
2159 }
2160
2161 if (color)
2162 {
2163 bool remove_primer = true;
2164 seg_str = convert_bp_to_color(org_seg_str, remove_primer);
2165 }
2166
2167 size_t to = 0;
2168 size_t seg_len = length(seg_str);
2169 size_t read_len = seg.support_read.size();
2170 if (read_len <= 0)
2171 to = seg_len - 2;
2172 else
2173 to = read_len - 2;
2174
2175 const size_t max_segment_len = 128;
2176 uint8_t left_mismatches[max_segment_len] = {0,};
2177 uint8_t right_mismatches[max_segment_len] = {0,};
2178
2179 if (max_segment_len < read_len)
2180 {
2181 fprintf(stderr, "Error: read len(%d) is greater than %d\n", (int)read_len, (int)max_segment_len);
2182 exit(-1);
2183 }
2184
2185 if (read_len == seg_len || seg.points_where == POINT_DIR_BOTH)
2186 {
2187 if (seg.points_where == POINT_DIR_RIGHT || seg.points_where == POINT_DIR_BOTH)
2188 {
2189 size_t num_mismatches = 0;
2190 for (size_t i = 0; i < read_len - 1; ++i)
2191 {
2192 if (seg_str[i] != seg.support_read[i])
2193 ++num_mismatches;
2194
2195 left_mismatches[i] = num_mismatches;
2196 if (num_mismatches > 2)
2197 {
2198 to = i;
2199 break;
2200 }
2201 }
2202 }
2203
2204 if (seg.points_where == POINT_DIR_LEFT || seg.points_where == POINT_DIR_BOTH)
2205 {
2206 size_t num_mismatches = 0;
2207 for (int i = read_len - 1; i >= 0; --i)
2208 {
2209 if (seg_str[i + (seg_len - read_len)] != seg.support_read[i])
2210 ++num_mismatches;
2211
2212 right_mismatches[i] = num_mismatches;
2213 if (num_mismatches > 2)
2214 break;
2215 }
2216 }
2217
2218 // daehwan
2219 if (bDebug)
2220 {
2221 cout << "antisense: " << (seg.antisense ? "-" : "+") << endl
2222 << seqan::infix(seg_str, 0, segment_length) << " "
2223 << seqan::infix(seg_str, length(seg_str) - segment_length, length(seg_str)) << endl
2224 << seg.support_read << endl
2225 << 0 << " - " << to << endl;
2226
2227 for (unsigned int i = 0; i < read_len; ++i)
2228 cout << (int)left_mismatches[i];
2229 cout << "\t";
2230
2231 for (unsigned int i = 0; i < read_len; ++i)
2232 cout << (int)right_mismatches[i];
2233 cout << endl;
2234 }
2235 }
2236
2237 if (seg.points_where == POINT_DIR_BOTH)
2238 {
2239 for (size_t i = 0; i <= to; ++i)
2240 {
2241 // Look at a slice of the reference without creating a copy.
2242 DnaString curr = seqan::infix(org_seg_str, i - left_color_offset, i + 2 - left_color_offset);
2243
2244 if ((!skip_fwd && curr == donor_dinuc) || (!skip_rev && curr == rev_acceptor_dinuc))
2245 {
2246 DnaString partner;
2247 if (curr == donor_dinuc)
2248 partner = acceptor_dinuc;
2249 else
2250 partner = rev_donor_dinuc;
2251
2252 uint8_t left_mismatch = 0;
2253 if (i > 0)
2254 left_mismatch = left_mismatches[i-1];
2255
2256 // daehwan
2257 if (bDebug)
2258 {
2259 cout << "i: " << i << endl
2260 << "mismatches: " << (int)left_mismatch
2261 << " - " << (int)right_mismatches[i] << endl;
2262 }
2263
2264 if (left_mismatch + right_mismatches[i] <= 2)
2265 {
2266 size_t pos = length(seg_str) - (read_len - i) - 2;
2267 if (partner == seqan::infix(org_seg_str, pos - left_color_offset, pos + 2 - left_color_offset))
2268 {
2269 if (curr == donor_dinuc)
2270 {
2271 motifs.fwd_donors.push_back(make_pair(seg.left + i, DnaSpliceStrings(0,0)));
2272 motifs.fwd_acceptors.push_back(make_pair(seg.left + pos, DnaSpliceStrings(0,0)));
2273 }
2274 else
2275 {
2276 motifs.rev_acceptors.push_back(make_pair(seg.left + i, DnaSpliceStrings(0,0)));
2277 motifs.rev_donors.push_back(make_pair(seg.left + pos, DnaSpliceStrings(0,0)));
2278 }
2279
2280 // daehwan
2281 if (bDebug)
2282 {
2283 cout << curr << ":" << partner << " added" << endl;
2284 }
2285 }
2286 }
2287 }
2288 }
2289 }
2290
2291 else if (seg.points_where == POINT_DIR_LEFT)
2292 {
2293 // A ref segment that "points left" is one that was flanked
2294 // on the right by a partial bowtie hit, indicating that we
2295 // should be looking for an intron to the left of the hit
2296 // In this seg, that means either an "AG" or an "AC"
2297 for (size_t i = 0; i <= to; ++i)
2298 {
2299 // Look at a slice of the reference without creating a copy.
2300 DnaString curr = seqan::infix(org_seg_str, i - left_color_offset, i + 2 - left_color_offset);
2301
2302 if (curr == acceptor_dinuc && !skip_fwd)
2303 motifs.fwd_acceptors.push_back(make_pair(seg.left + i, DnaSpliceStrings(0,0)));
2304 else if (curr == rev_donor_dinuc && !skip_rev)
2305 motifs.rev_donors.push_back(make_pair(seg.left + i, DnaSpliceStrings(0,0)));
2306 }
2307 }
2308 else
2309 {
2310 // A right pointing ref seg wants either a "GT" or a "CT"
2311 for (size_t i = 0; i <= to; ++i)
2312 {
2313 // Look at a slice of the reference without creating a copy.
2314 DnaString curr = seqan::infix(org_seg_str, i - left_color_offset, i + 2 - left_color_offset);
2315
2316 if (curr == donor_dinuc && !skip_fwd)
2317 motifs.fwd_donors.push_back(make_pair(seg.left + i, DnaSpliceStrings(0,0)));
2318 else if (curr == rev_acceptor_dinuc && !skip_rev)
2319 motifs.rev_acceptors.push_back(make_pair(seg.left + i, DnaSpliceStrings(0,0)));
2320 }
2321 }
2322 }
2323
2324 if (talkative)
2325 {
2326 fprintf(stderr, "reporting synthetic splice junctions...\n");
2327 }
2328 for (MotifMap::iterator motif_itr = ims.begin(); motif_itr != ims.end(); ++motif_itr)
2329 {
2330 uint32_t ref_id = motif_itr->first;
2331
2332 RefSequenceTable::Sequence* ref_str = rt.get_seq(ref_id);
2333
2334 if (!ref_str)
2335 err_die("Error: couldn't get ref string for %u\n", ref_id);
2336
2337 if (talkative)
2338 fprintf(stderr, "Examining donor-acceptor pairings in %s\n", rt.get_name(ref_id));
2339 IntronMotifs& motifs = motif_itr->second;
2340
2341 if (!all_both)
2342 motifs.unique();
2343
2344 //motifs.attach_mer_counts(*ref_str);
2345 motifs.attach_mers(*ref_str);
2346
2347 vector<pair<size_t, DnaSpliceStrings> >& fwd_donors = motifs.fwd_donors;
2348 vector<pair<size_t, DnaSpliceStrings> >& fwd_acceptors = motifs.fwd_acceptors;
2349 vector<pair<size_t, DnaSpliceStrings> >& rev_acceptors = motifs.rev_acceptors;
2350 vector<pair<size_t, DnaSpliceStrings> >& rev_donors = motifs.rev_donors;
2351
2352 //const char* ref_name = rt.get_name(motif_itr->second.ref_id);
2353
2354 JunctionRecorder recorder;
2355 recorder.record(ref_id,
2356 fwd_donors,
2357 fwd_acceptors,
2358 false,
2359 juncs,
2360 min_intron,
2361 max_intron,
2362 max_juncs,
2363 half_splice_mer_len);
2364
2365 recorder.record(ref_id,
2366 rev_acceptors,
2367 rev_donors,
2368 true,
2369 juncs,
2370 min_intron,
2371 max_intron,
2372 max_juncs,
2373 half_splice_mer_len);
2374 }
2375 //fprintf(stderr, "Found %d total splices\n", num_juncs);
2376 }
2377
2378
2379 /**
2380 * Performs a simple global alignment.
2381 * This function will perform a restricted global alignment. The restriction is that only one insertion/deletion
2382 * is allowed in the final alignment.
2383 * @param shortSequence The short sequence to be aligned.
2384 * @param leftReference The left end of the reference to be aligned, must be exactly as long as the short sequence
2385 * @param leftReference The right end of the reference to be aligned, must be exactly as long as the short sequenc
2386 * @param insertPosition This will contain the 0-based index of the first position in the shorter sequence after the insertion/deletion. A value of -1 indicates that the alignment could not be performed.
2387 * @param mismatchCount This will contain the number of mismatches in the optimal restricted global alignment. The number and length of insertions/deletions is fixed. A value of -1 indicates that the alignment could not be performed.
2388 */
2389 void simpleSplitAlignment(seqan::String<char>& shorterSequence,
2390 seqan::String<char>& leftReference,
2391 seqan::String<char>& rightReference,
2392 int& insertPosition,
2393 int& mismatchCount)
2394 {
2395 /*
2396 * In this restricted alignment, we already know the length and number (1) of insertions/deletions.
2397 * We simply need to know where to put it. Do a linear scan through sequence counting the number of induced
2398 * errors before and after putting the insertion at each sequence.
2399 */
2400
2401 /*
2402 * Note that we could have a case, where both the alignment and the read have the unknonw
2403 * nucleotide ('N') and we don't want to reward cases where these characters match
2404 */
2405 vector<unsigned short> beforeErrors(seqan::length(shorterSequence));
2406 for(int idx = seqan::length(shorterSequence) - 1; idx >= 0; idx -= 1){
2407 unsigned short prevCount = 0;
2408 /*
2409 * We guarentee idx >= 0, so cast to hide the compiler
2410 * warning here
2411 */
2412 if(((size_t)idx) < seqan::length(shorterSequence) - 1){
2413 prevCount = beforeErrors.at(idx + 1);
2414 }
2415 unsigned short currentMismatch = 0;
2416 if(rightReference[idx] == 'N' || shorterSequence[idx] == 'N' || rightReference[idx] != shorterSequence[idx]){
2417 currentMismatch = 1;
2418 }
2419 beforeErrors.at(idx) = prevCount + currentMismatch;
2420 }
2421
2422
2423 vector<unsigned short> afterErrors(seqan::length(shorterSequence));
2424 for(size_t idx = 0; idx < seqan::length(shorterSequence) ; idx += 1){
2425 unsigned short prevCount = 0;
2426 if(idx > 0){
2427 prevCount = afterErrors.at(idx - 1);
2428 }
2429 unsigned short currentMismatch = 0;
2430 if(leftReference[idx] == 'N' || shorterSequence[idx] == 'N' || leftReference[idx] != shorterSequence[idx]){
2431 currentMismatch = 1;
2432 }
2433 afterErrors.at(idx) = prevCount + currentMismatch;
2434 }
2435
2436
2437 mismatchCount = seqan::length(shorterSequence) + 1;
2438 insertPosition = -1;
2439
2440 /*
2441 * Technically, we could allow the insert position to be at the end or beginning of the sequence,
2442 * but we are disallowing it here
2443 */
2444 for(size_t currentInsertPosition = 1; currentInsertPosition < seqan::length(shorterSequence); currentInsertPosition += 1){
2445 size_t errorCount = beforeErrors.at(currentInsertPosition) + afterErrors.at(currentInsertPosition - 1);
2446 if(((int)errorCount) < mismatchCount){
2447 mismatchCount = (int)errorCount;
2448 insertPosition = currentInsertPosition;
2449 }
2450 }
2451 return;
2452 }
2453
2454
2455 /**
2456 * Try to detect a small insertion.
2457 * This code will try to identify a small insertion based on the ungapped alignment of two neighboring
2458 * segments. The general idea is to try to realign the local region, and see if we can reduce the
2459 * number of errors. Note that the function makes use of the global parameter "max_insertion_length" to limit the maximum
2460 * size of a detected insertion.
2461 * @param rt Sequence table used to lookup sequence information
2462 * @param leftHit The alignment of the left segment. Note that the leftHit must have a left position less than that of the right hit.
2463 * @param rightHit The alignment of the right segment. Note that the rightHit must have a left position greater than that of the left hit.
2464 * @param insertions If an insertion is sucessfully detected, it will be added to this set
2465 */
2466 void detect_small_insertion(RefSequenceTable& rt,
2467 seqan::String<char>& read_sequence,
2468 BowtieHit& leftHit,
2469 BowtieHit& rightHit,
2470 std::set<Insertion>& insertions)
2471 {
2472
2473 RefSequenceTable::Sequence* ref_str = rt.get_seq(leftHit.ref_id());
2474 if(!ref_str){
2475 fprintf(stderr, "Error in accessing sequence record\n");
2476 }else{
2477 size_t read_length = seqan::length(read_sequence);
2478 int begin_offset = 0;
2479 int end_offset = 0;
2480
2481 if(color){
2482 if(leftHit.antisense_align())
2483 end_offset = 1;
2484 else
2485 begin_offset = -1;
2486 }
2487
2488 if(leftHit.left() + begin_offset < 0)
2489 return;
2490
2491 /*
2492 * If there is in fact a deletion, we are expecting the genomic sequence to be shorter than
2493 * the actual read sequence
2494 */
2495 int discrepancy = read_length - (rightHit.right() - leftHit.left());
2496 DnaString genomic_sequence_temp = seqan::infix(*ref_str, leftHit.left() + begin_offset, rightHit.right() + end_offset);
2497 String<char> genomic_sequence;
2498 assign(genomic_sequence, genomic_sequence_temp);
2499
2500 if(color)
2501 genomic_sequence = convert_bp_to_color(genomic_sequence, true);
2502
2503 String<char> left_read_sequence = seqan::infix(read_sequence, 0, 0 + seqan::length(genomic_sequence));
2504 String<char> right_read_sequence = seqan::infix(read_sequence, read_length - seqan::length(genomic_sequence), read_length);
2505
2506 int bestInsertPosition = -1;
2507 int minErrors = -1;
2508 simpleSplitAlignment(genomic_sequence, left_read_sequence, right_read_sequence, bestInsertPosition, minErrors);
2509
2510 /*
2511 * Need to decide if the insertion is suitably improves the alignment
2512 */
2513 /*
2514 * If these two segments anchors constitue the entire read, then we require
2515 * that this alignment actually improve the number of errors observed in the alignment
2516 * Otherwise, it is OK as long as the number of errors doesn't increase.
2517 */
2518 int adjustment = 0;
2519 if(leftHit.read_len() + rightHit.read_len() >= (int)read_length){
2520 adjustment = -1;
2521 }
2522 if(minErrors <= (leftHit.edit_dist()+rightHit.edit_dist()+adjustment)){
2523 String<char> insertedSequence = seqan::infix(left_read_sequence, bestInsertPosition, bestInsertPosition + discrepancy);
2524 if(color)
2525 insertedSequence = convert_color_to_bp(genomic_sequence_temp[bestInsertPosition - begin_offset + end_offset - 1], insertedSequence);
2526
2527 insertions.insert(Insertion(leftHit.ref_id(),
2528 leftHit.left() + bestInsertPosition - 1 + end_offset,
2529 seqan::toCString(insertedSequence)));
2530 }
2531 }
2532 return;
2533 }
2534
2535 /**
2536 * Try to detect a small deletion.
2537 * This code will try to identify a small deletion based on the ungapped alignment of two neighboring
2538 * segments. The general idea is to try to realign the local region, and see if we can reduce the
2539 * number of errors. Note that the function makes use of the global parameter "max_deletion_length" to limit the maximum
2540 * size of a detected deletion.
2541 * @param rt Sequence table used to lookup sequence information
2542 * @param leftHit The alignment of the left segment. Note that the leftHit must have a left position less than that of the right hit.
2543 * @param rightHit The alignment of the right segment. Note that the rightHit must have a left position greater than that of the left hit.
2544 * @param deletion_juncs If a deletion is sucessfully detected, it will be added to this set
2545 */
2546 void detect_small_deletion(RefSequenceTable& rt,
2547 seqan::String<char>& read_sequence,
2548 BowtieHit& leftHit,
2549 BowtieHit& rightHit,
2550 std::set<Deletion>& deletions)
2551 {
2552
2553 RefSequenceTable::Sequence* ref_str = rt.get_seq(leftHit.ref_id());
2554 if(!ref_str){
2555 fprintf(stderr, "Error in accessing sequence record\n");
2556 }else{
2557 int begin_offset = 0;
2558 int end_offset = 0;
2559
2560 if(color){
2561 if(leftHit.antisense_align())
2562 end_offset = 1;
2563 else
2564 begin_offset = -1;
2565 }
2566
2567 if(leftHit.left() + begin_offset < 0)
2568 return;
2569
2570 size_t read_length = seqan::length(read_sequence);
2571 if(rightHit.right() + read_length + begin_offset < 0 )
2572 return;
2573
2574 int discrepancy = (rightHit.right() - leftHit.left()) - read_length;
2575 Dna5String leftGenomicSequence_temp = seqan::infix(*ref_str, leftHit.left() + begin_offset, leftHit.left() + read_length + end_offset);
2576 Dna5String rightGenomicSequence_temp = seqan::infix(*ref_str, rightHit.right() - read_length + begin_offset, rightHit.right() + end_offset);
2577
2578 if (length(leftGenomicSequence_temp) < read_length || length(rightGenomicSequence_temp) < read_length)
2579 return;
2580
2581 String<char> leftGenomicSequence;
2582 assign(leftGenomicSequence, leftGenomicSequence_temp);
2583
2584 String<char> rightGenomicSequence;
2585 assign(rightGenomicSequence, rightGenomicSequence_temp);
2586
2587 if(color){
2588 leftGenomicSequence = convert_bp_to_color(leftGenomicSequence, true);
2589 rightGenomicSequence = convert_bp_to_color(rightGenomicSequence, true);
2590 }
2591
2592 int bestInsertPosition = -1;
2593 int minErrors = -1;
2594
2595 simpleSplitAlignment(read_sequence, leftGenomicSequence, rightGenomicSequence, bestInsertPosition, minErrors);
2596
2597 /*
2598 * Need to decide if the deletion is suitably improves the alignment
2599 */
2600 int adjustment = 0;
2601
2602 /*
2603 * If these two segments anchors constitue the entire read, then we require
2604 * that this alignment actually improve the number of errors observed in the alignment
2605 * Otherwise, it is OK as long as the number of errors doesn't increase.
2606 */
2607 if(leftHit.read_len() + rightHit.read_len() >= (int)read_length){
2608 adjustment = -1;
2609 }
2610 if(minErrors <= (leftHit.edit_dist()+rightHit.edit_dist()+adjustment)){
2611 deletions.insert(Deletion(leftHit.ref_id(),
2612 leftHit.left() + bestInsertPosition - 1 + end_offset,
2613 leftHit.left() + bestInsertPosition + discrepancy + end_offset,
2614 false));
2615 }
2616 }
2617 return;
2618 }
2619
2620
2621 void find_insertions_and_deletions(RefSequenceTable& rt,
2622 FILE* reads_file,
2623 vector<HitsForRead>& hits_for_read,
2624 std::set<Deletion>& deletions,
2625 std::set<Insertion>& insertions){
2626
2627 if(hits_for_read.empty()){
2628 return;
2629 }
2630
2631 size_t last_segment = hits_for_read.size()-1;
2632 size_t first_segment = 0;
2633 if(last_segment == first_segment){
2634 return;
2635 }
2636
2637 /*
2638 * We can check up front whether the first or last element is empty
2639 * and avoid doing any more work. Note that the following code requires
2640 * that there be at least one elment in each
2641 */
2642 if(hits_for_read[first_segment].hits.empty() || hits_for_read[last_segment].hits.empty()){
2643 return;
2644 }
2645
2646 char read_name[1024];
2647 char read_seq[1024];
2648 char read_alt_name[1024];
2649 char read_quals[1024];
2650
2651 /*
2652 * Need to identify the appropriate insert id for this group of reads
2653 */
2654
2655 bool got_read = get_read_from_stream(hits_for_read.back().insert_id,
2656 reads_file,
2657 FASTQ,
2658 false,
2659 read_name,
2660 read_seq,
2661 read_alt_name,
2662 read_quals);
2663 if(!got_read){
2664 return;
2665 }
2666
2667 /*
2668 * Work through all combinations of mappings for the first and last segment to see if any are indicative
2669 * of a small insertions or deletion
2670 */
2671 HitsForRead& left_segment_hits = hits_for_read[first_segment];
2672 HitsForRead& right_segment_hits = hits_for_read[last_segment];
2673
2674 /*
2675 * If either of the segment match lists is empty, we could try
2676 * to be smarter and work our way in until we find good a segment
2677 * match; however, won't do that for noe.
2678 */
2679 if(left_segment_hits.hits.empty() || right_segment_hits.hits.empty()){
2680 return;
2681 }
2682
2683 seqan::String<char> fullRead, rcRead;
2684 if(color){
2685 fullRead = read_seq + 1;
2686 rcRead = fullRead;
2687 seqan::reverseInPlace(rcRead);
2688 }else{
2689 fullRead = read_seq;
2690 rcRead = read_seq;
2691 seqan::convertInPlace(rcRead, seqan::FunctorComplement<Dna>());
2692 seqan::reverseInPlace(rcRead);
2693 }
2694
2695 size_t read_length = seqan::length(fullRead);
2696 for(size_t left_segment_index = 0; left_segment_index < left_segment_hits.hits.size(); left_segment_index++){
2697 for(size_t right_segment_index = 0; right_segment_index < right_segment_hits.hits.size(); right_segment_index++){
2698 BowtieHit* leftHit = &left_segment_hits.hits[left_segment_index];
2699 BowtieHit* rightHit = &right_segment_hits.hits[right_segment_index];
2700 /*
2701 * Now we have found a pair of segment hits to investigate. Need to ensure
2702 * that
2703 * 1. the alignment orientation is consistent
2704 * 2. the distance separation is in the appropriate range
2705 * 3. Both hits are aligned to the same contig
2706 */
2707 if(leftHit->ref_id() != rightHit->ref_id()){
2708 continue;
2709 }
2710
2711 if(leftHit->antisense_align() != rightHit->antisense_align()){
2712 continue;
2713 }
2714
2715 seqan::String<char>* modifiedRead = &fullRead;
2716 /*
2717 * If we are dealing with an antisense alignment, then the left
2718 * read will actually be on the right, fix this now, to simplify
2719 * the rest of the logic, in addition, we will need to use the reverse
2720 * complement of the read sequence
2721 */
2722 if(leftHit->antisense_align()){
2723 BowtieHit * tmp = leftHit;
2724 leftHit = rightHit;
2725 rightHit = tmp;
2726 modifiedRead = &rcRead;
2727 }
2728
2729 size_t apparent_length = rightHit->right() - leftHit->left();
2730 int length_discrepancy = apparent_length - read_length;
2731 if(length_discrepancy > 0 && length_discrepancy <= (int)max_deletion_length){
2732 /*
2733 * Search for a deletion
2734 */
2735 detect_small_deletion(rt, *modifiedRead, *leftHit, *rightHit, deletions);
2736 }
2737 if(length_discrepancy < 0 && length_discrepancy >= -(int)max_insertion_length){
2738 /*
2739 * Search for an insertion
2740 */
2741 detect_small_insertion(rt, *modifiedRead, *leftHit, *rightHit, insertions);
2742 }
2743 }
2744 }
2745 }
2746
2747
2748 void find_gaps(RefSequenceTable& rt,
2749 FILE* reads_file,
2750 vector<HitsForRead>& hits_for_read,
2751 std::set<Junction, skip_count_lt>& seg_juncs,
2752 eREAD read)
2753 {
2754 if (hits_for_read.empty())
2755 return;
2756
2757 size_t last_non_empty = hits_for_read.size() - 1;
2758 while(last_non_empty >= 0 && hits_for_read[last_non_empty].hits.empty())
2759 --last_non_empty;
2760
2761 hits_for_read.resize(last_non_empty + 1);
2762
2763 char read_name[1024];
2764 char read_seq[1024];
2765 char read_alt_name[1024];
2766 char read_quals[1024];
2767 bool got_read = get_read_from_stream(hits_for_read[last_non_empty].insert_id,
2768 reads_file,
2769 FASTQ,
2770 false,
2771 read_name,
2772 read_seq,
2773 read_alt_name,
2774 read_quals);
2775
2776 if (!got_read)
2777 return;
2778
2779 vector<RefSeg> expected_don_acc_windows;
2780 string seq(read_seq);
2781
2782 for (size_t s = 0; s < hits_for_read.size(); ++s)
2783 {
2784 HitsForRead& curr = hits_for_read[s];
2785 for (size_t h = 0; h < curr.hits.size(); ++h)
2786 {
2787 bool found_right_partner = s == hits_for_read.size() - 1;
2788 bool found_distant_right_partner = false;
2789 bool found_right_right_partner = false;
2790 BowtieHit& bh = curr.hits[h];
2791
2792 // "dr" is distant right partner
2793 // "rr" is right of right partner
2794 BowtieHit dr_bh, rr_bh;
2795 if (s < hits_for_read.size() - 1)
2796 {
2797 // Look for a right partner for the current hit
2798 HitsForRead& right = hits_for_read[s + 1];
2799
2800 for (size_t r = 0; r < right.hits.size(); ++r)
2801 {
2802 BowtieHit& rh = right.hits[r];
2803 if (bh.antisense_align() != rh.antisense_align() || bh.ref_id() != rh.ref_id())
2804 continue;
2805
2806 if ((bh.antisense_align() && rh.right() == bh.left()) ||
2807 (!bh.antisense_align() && bh.right() == rh.left() ))
2808 {
2809 found_right_partner = true;
2810 break;
2811 }
2812
2813 int dist = 0;
2814 if (bh.antisense_align())
2815 dist = bh.left() - rh.right();
2816 else
2817 dist = rh.left() - bh.right();
2818
2819 if (dist >= min_segment_intron_length && dist < (int)max_segment_intron_length)
2820 {
2821 found_distant_right_partner = true;
2822 dr_bh = rh;
2823 }
2824 }
2825 }
2826
2827 if (!found_right_partner && s < hits_for_read.size() - 2)
2828 {
2829 // Look for a right of right partner for the current hit
2830 HitsForRead& right_right = hits_for_read[s + 2];
2831
2832 for (size_t r = 0; r < right_right.hits.size(); ++r)
2833 {
2834 BowtieHit& rrh = right_right.hits[r];
2835 if (bh.antisense_align() != rrh.antisense_align() || bh.ref_id() != rrh.ref_id())
2836 continue;
2837
2838 int dist = 0;
2839 if (bh.antisense_align())
2840 dist = bh.left() - rrh.right();
2841 else
2842 dist = rrh.left() - bh.right();
2843
2844 if (dist >= min_segment_intron_length + segment_length && dist < (int)max_segment_intron_length + segment_length)
2845 {
2846 found_right_right_partner = true;
2847 rr_bh = rrh;
2848 break;
2849 }
2850 }
2851 }
2852
2853 if (!found_right_partner && (found_distant_right_partner || found_right_right_partner))
2854 {
2855 size_t color_offset = color ? 1 : 0;
2856
2857 string support_read;
2858 if (found_distant_right_partner)
2859 support_read = seq.substr(color_offset + (s+1) * segment_length - 5, 10);
2860 else
2861 support_read = seq.substr(color_offset + (s+1) * segment_length - 5, segment_length + 10);
2862
2863 BowtieHit& d_bh = found_distant_right_partner ? dr_bh : rr_bh;
2864 if (!bh.antisense_align())
2865 {
2866 RefSeg right_seg(bh.ref_id(), POINT_DIR_BOTH, bh.antisense_align(), read, 0, 0, support_read);
2867 right_seg.left = max(0, bh.right() - 5);
2868 right_seg.right = d_bh.left() + 5;
2869 expected_don_acc_windows.push_back(right_seg);
2870 }
2871 else
2872 {
2873 if (color)
2874 reverse(support_read.begin(), support_read.end());
2875 else
2876 reverse_complement(support_read);
2877
2878 RefSeg left_seg(bh.ref_id(), POINT_DIR_BOTH, bh.antisense_align(), read, 0, 0, support_read);
2879 left_seg.left = d_bh.right() - 5;
2880 left_seg.right = bh.left() + 5; // num allowed bowtie mismatches
2881 expected_don_acc_windows.push_back(left_seg);
2882 }
2883
2884 // daehwan
2885 if (bDebug)
2886 {
2887 cout << (bh.antisense_align() ? "-" : "+") << endl
2888 << seq << endl
2889 << "(" << s << ") - " << support_read << endl;
2890 }
2891 }
2892 }
2893 }
2894
2895 juncs_from_ref_segs<RecordSegmentJuncs>(rt,
2896 expected_don_acc_windows,
2897 seg_juncs,
2898 "GT",
2899 "AG",
2900 max_segment_intron_length,
2901 min_segment_intron_length,
2902 max_seg_juncs,
2903 false,
2904 0);
2905
2906 juncs_from_ref_segs<RecordSegmentJuncs>(rt,
2907 expected_don_acc_windows,
2908 seg_juncs,
2909 "GC",
2910 "AG",
2911 max_segment_intron_length,
2912 min_segment_intron_length,
2913 max_seg_juncs,
2914 false,
2915 0);
2916
2917 juncs_from_ref_segs<RecordSegmentJuncs>(rt,
2918 expected_don_acc_windows,
2919 seg_juncs,
2920 "AT",
2921 "AC",
2922 max_segment_intron_length,
2923 min_segment_intron_length,
2924 max_seg_juncs,
2925 false,
2926 0);
2927 }
2928
2929
2930 MerTable mer_table;
2931
2932 int seed_alignments = 0;
2933 int microaligned_segs = 0;
2934
2935 map<RefSeg, vector<string>* > microexon_windows;
2936
2937 bool overlap_in_genome(int ll, int lr, int rl, int rr)
2938 {
2939 if (ll >= rl && ll < rr)
2940 return true;
2941 if (lr > rl && lr < rr)
2942 return true;
2943 if (rl >= ll && rl < lr)
2944 return true;
2945 if (rr > ll && rr < lr)
2946 return true;
2947 return false;
2948 }
2949
2950 void add_to_microexon_windows(uint32_t ref_id,
2951 int left_boundary,
2952 int right_boundary,
2953 const string& dna_str,
2954 eREAD read)
2955 {
2956 RefSeg left_dummy(ref_id, POINT_DIR_DONTCARE, false, read, left_boundary, right_boundary);
2957 RefSeg right_dummy(ref_id, POINT_DIR_DONTCARE, false, read, right_boundary, right_boundary + 1);
2958
2959 map<RefSeg, vector<string>* >::iterator lb = microexon_windows.lower_bound(left_dummy);
2960 map<RefSeg, vector<string>* >::iterator ub = microexon_windows.lower_bound(right_dummy);
2961 vector<string>* new_vec = NULL;
2962 if (lb == microexon_windows.end())
2963 {
2964 microexon_windows.insert(make_pair(left_dummy, new vector<string>(1, dna_str)));
2965 return;
2966 }
2967
2968 map<RefSeg, vector<string>* >::iterator first_to_be_erased = microexon_windows.end();
2969 map<RefSeg, vector<string>* >::iterator last_to_be_erased = ub;
2970
2971 while (lb != ub)
2972 {
2973 // everyone in this range that overlaps with the new interval needs to
2974 // be merged together.
2975 if (overlap_in_genome(lb->first.left, lb->first.right, left_boundary, right_boundary))
2976 {
2977 if (!new_vec)
2978 new_vec = new vector<string>();
2979 if (first_to_be_erased == microexon_windows.end())
2980 first_to_be_erased = lb;
2981 left_dummy.left = min(lb->first.left, left_boundary);
2982 left_dummy.right = max(lb->first.right, right_boundary);
2983
2984 new_vec->insert(new_vec->end(), lb->second->begin(), lb->second->end());
2985 delete lb->second;
2986 }
2987 else if (first_to_be_erased != microexon_windows.end())
2988 {
2989 last_to_be_erased = lb;
2990 }
2991
2992 ++lb;
2993 }
2994
2995 if (first_to_be_erased != microexon_windows.end())
2996 {
2997 microexon_windows.erase(first_to_be_erased, last_to_be_erased);
2998 }
2999
3000 if (!new_vec)
3001 {
3002 // never found an overlapping window, so just add this one and bail
3003 microexon_windows.insert(make_pair(left_dummy, new vector<string>(1, dna_str)));
3004 return;
3005 }
3006 else
3007 {
3008 new_vec->push_back(dna_str);
3009 microexon_windows.insert(make_pair(left_dummy, new_vec));
3010 return;
3011 }
3012
3013 }
3014
3015 void align_microexon_segs(RefSequenceTable& rt,
3016 std::set<Junction, skip_count_lt>& juncs,
3017 int max_juncs,
3018 int half_splice_mer_len)
3019 {
3020 int num_segments = 0;
3021 for (map<RefSeg, vector<string>* >::iterator itr = microexon_windows.begin();
3022 itr != microexon_windows.end(); ++itr)
3023 {
3024 vector<string>& unaligned_segments = *itr->second;
3025 num_segments += unaligned_segments.size();
3026 }
3027
3028 fprintf(stderr, "Aligning %d microexon segments in %lu windows\n",
3029 num_segments, (long unsigned int)microexon_windows.size());
3030
3031 extensions.clear();
3032
3033 size_t splice_mer_len = 2 * half_splice_mer_len;
3034 size_t mer_table_size = 1 << ((splice_mer_len)<<1);
3035
3036 extensions.resize(mer_table_size);
3037
3038 int window_num = 0;
3039 for (map<RefSeg, vector<string>* >::iterator itr = microexon_windows.begin();
3040 itr != microexon_windows.end(); ++itr)
3041 {
3042 window_num++;
3043 if ((window_num % 100) == 0)
3044 fprintf(stderr, "\twindow %d\n",window_num);
3045
3046
3047 stringstream ss(stringstream::in | stringstream::out);
3048
3049 for (size_t j = 0; j < extensions.size(); ++j)
3050 {
3051 extensions[j].clear();
3052 }
3053
3054 vector<string>& unaligned_segments = *itr->second;
3055
3056 for (size_t j = 0; j < unaligned_segments.size(); ++j)
3057 {
3058 stringstream ss(stringstream::in | stringstream::out);
3059 string s;
3060 //cerr << w.unaligned_segments[j];
3061 ss << unaligned_segments[j];
3062 ss >> s;
3063
3064 store_read_extensions(extensions,
3065 half_splice_mer_len,
3066 half_splice_mer_len,
3067 s,
3068 false);
3069 }
3070
3071 vector<RefSeg> segs;
3072 segs.push_back(itr->first);
3073 RefSeg r = itr->first;
3074 r.points_where = POINT_DIR_LEFT;
3075 segs.push_back(r);
3076
3077 juncs_from_ref_segs<RecordExtendableJuncs>(rt,
3078 segs,
3079 juncs,
3080 "GT",
3081 "AG",
3082 max_microexon_stretch,
3083 min_coverage_intron_length,
3084 max_juncs,
3085 false,
3086 half_splice_mer_len);
3087 num_segments += unaligned_segments.size();
3088 delete itr->second;
3089 }
3090 fprintf(stderr, "Checked %d segments against %lu windows for microexon junctions\n",
3091 num_segments, (long unsigned int)microexon_windows.size());
3092 fprintf(stderr, "Found %ld potential microexon junctions\n", (long int)juncs.size());
3093 }
3094
3095 /*
3096 * Easy guys ... this function puts its pants on just like the rest of you -- one leg
3097 * at a time. Except, once its pants are on, it makes gold records. Alright, here we
3098 * go.
3099 */
3100
3101 void look_for_hit_group(RefSequenceTable& rt,
3102 //FILE* reads_file,
3103 FZPipe& reads_file,
3104 FZPipe& reads_file_for_segment_search,
3105 FZPipe& reads_file_for_indel_discovery,
3106 ReadTable& unmapped_reads,
3107 vector<HitStream>& seg_files,
3108 int curr_file,
3109 uint64_t insert_id,
3110 vector<HitsForRead>& hits_for_read,
3111 std::set<Junction, skip_count_lt>& juncs,
3112 std::set<Deletion>& deletions,
3113 std::set<Insertion>& insertions,
3114 eREAD read)
3115
3116 {
3117 HitStream& hit_stream = seg_files[curr_file];
3118 HitsForRead hit_group;
3119 uint32_t order = unmapped_reads.observation_order(insert_id);
3120 int seq_key_len = min((int)min_anchor_len, 6);
3121 while(true) {
3122 uint64_t next_group_id = hit_stream.next_group_id();
3123 uint32_t next_order = unmapped_reads.observation_order(next_group_id);
3124 // If we would have seen the hits by now, stop looking in this stream,
3125 // but forward the search if possible.
3126 if (order < next_order || next_group_id == 0) {
3127 if (curr_file > 0) {
3128 look_for_hit_group(rt,
3129 reads_file,
3130 reads_file_for_segment_search,
3131 reads_file_for_indel_discovery,
3132 unmapped_reads,
3133 seg_files,
3134 curr_file - 1,
3135 insert_id,
3136 hits_for_read,
3137 juncs,
3138 deletions,
3139 insertions,
3140 read);
3141 }
3142 else
3143 if (!no_microexon_search) {
3144 //microexon search
3145 char read_name[1024];
3146 char read_seq[1024];
3147 char read_alt_name[1024];
3148 char read_quals[1024];
3149 // The hits are missing for the leftmost segment, which means
3150 // we should try looking for junctions via seed and extend
3151 // using it (helps find junctions to microexons).
3152 bool got_read = get_read_from_stream(insert_id,
3153 reads_file.file,
3154 FASTQ,
3155 false,
3156 read_name,
3157 read_seq,
3158 read_alt_name,
3159 read_quals);
3160 if (!got_read) {
3161 fprintf(stderr, "Warning: could not get read with insert_id=%d\n", (int)insert_id);
3162 break; //exit loop
3163 }
3164 string fwd_read = read_seq;
3165 if (color) // remove the primer and the adjacent color
3166 fwd_read.erase(0, 2);
3167 // make sure there are hits for all the other segs, all the
3168 // way to the root (right-most) one.
3169 int empty_seg = 0;
3170 for (size_t h = 1; h < hits_for_read.size(); ++h) {
3171 if (hits_for_read[h].hits.empty())
3172 empty_seg = h;
3173 }
3174 // if not, no microexon alignment for this segment
3175 if (empty_seg != 0)
3176 break;
3177 fwd_read = fwd_read.substr(0, segment_length);
3178 string rev_read = fwd_read;
3179 //check the reverse
3180 if (color) reverse(rev_read.begin(), rev_read.end());
3181 else reverse_complement(rev_read);
3182 for (size_t h = 0; h < hits_for_read[empty_seg + 1].hits.size(); ++h) {
3183 const BowtieHit& bh = hits_for_read[empty_seg + 1].hits[h];
3184 RefSequenceTable::Sequence* ref_str = rt.get_seq(bh.ref_id());
3185 if (ref_str == NULL)
3186 continue;
3187 int ref_len = length(*ref_str);
3188 int left_boundary;
3189 int right_boundary;
3190 bool antisense = bh.antisense_align();
3191 vector<BowtieHit> empty_seg_hits;
3192 seed_alignments++;
3193 if (antisense) {
3194 left_boundary = max(0, bh.right() - (int)min_anchor_len);
3195 right_boundary = min(ref_len - 2,
3196 left_boundary + max_microexon_stretch);
3197 if (right_boundary - left_boundary < 2 * seq_key_len)
3198 continue;
3199 microaligned_segs++;
3200 add_to_microexon_windows(bh.ref_id(), left_boundary, right_boundary, rev_read, read);
3201 }
3202 else {
3203 right_boundary = min(ref_len - 2,
3204 bh.left() + (int)min_anchor_len);
3205 left_boundary = max(0, right_boundary - max_microexon_stretch);
3206 if (right_boundary - left_boundary < 2 * seq_key_len)
3207 continue;
3208 microaligned_segs++;
3209 add_to_microexon_windows(bh.ref_id(), left_boundary, right_boundary, fwd_read, read);
3210 }
3211 } //for h
3212 } // !no_microexon_search
3213 break;
3214 }
3215 else if (hit_stream.next_read_hits(hit_group)) {
3216 // if we found hits for the target group in the left stream,
3217 // add them to the accumulating vector and continue the search
3218 if (hit_group.insert_id == insert_id) {
3219 hits_for_read[curr_file] = hit_group;
3220 if (curr_file > 0)
3221 // we need to look left (recursively) for the group we
3222 // just read for this stream.
3223 look_for_hit_group(rt,
3224 reads_file,
3225 reads_file_for_segment_search,
3226 reads_file_for_indel_discovery,
3227 unmapped_reads,
3228 seg_files,
3229 curr_file - 1,
3230 insert_id,
3231 hits_for_read,
3232 juncs,
3233 deletions,
3234 insertions,
3235 read);
3236 break;
3237 } //same insert_id (group)
3238 else if (curr_file > 0) {
3239 // different group, we need to start a whole new
3240 // search for it, with a whole new vector of hits.
3241 vector<HitsForRead> hits_for_new_read(seg_files.size());
3242 hits_for_new_read[curr_file] = hit_group;
3243 look_for_hit_group(rt,
3244 reads_file,
3245 reads_file_for_segment_search,
3246 reads_file_for_indel_discovery,
3247 unmapped_reads,
3248 seg_files,
3249 curr_file - 1,
3250 hit_group.insert_id,
3251 hits_for_new_read,
3252 juncs,
3253 deletions,
3254 insertions,
3255 read);
3256 find_insertions_and_deletions(rt, reads_file_for_indel_discovery.file, hits_for_new_read, deletions, insertions);
3257 find_gaps(rt, reads_file_for_segment_search.file, hits_for_new_read, juncs, read);
3258 } //different group
3259 }//got next group
3260 } //while loop
3261 }
3262
3263
3264
3265 bool process_next_hit_group(RefSequenceTable& rt,
3266 //FILE* reads_file,
3267 FZPipe& reads_file,
3268 FZPipe& reads_file_for_segment_search,
3269 FZPipe& reads_file_for_indel_discovery,
3270 ReadTable& unmapped_reads,
3271 vector<HitStream>& seg_files,
3272 size_t curr_file,
3273 std::set<Junction, skip_count_lt>& juncs,
3274 std::set<Deletion>& deletions,
3275 std::set<Insertion>& insertions,
3276 eREAD read)
3277 {
3278 HitStream& curr = seg_files[curr_file];
3279 HitsForRead hit_group;
3280 bool result = curr.next_read_hits(hit_group);
3281 vector<HitsForRead> hits_for_read(seg_files.size());
3282 hits_for_read.back() = hit_group;
3283
3284 look_for_hit_group(rt,
3285 reads_file,
3286 reads_file_for_segment_search,
3287 reads_file_for_indel_discovery,
3288 unmapped_reads,
3289 seg_files,
3290 (int)curr_file - 1,
3291 hit_group.insert_id,
3292 hits_for_read,
3293 juncs,
3294 deletions,
3295 insertions,
3296 read);
3297
3298 if (result) {
3299 find_insertions_and_deletions(rt, reads_file_for_indel_discovery.file, hits_for_read, deletions, insertions);
3300 find_gaps(rt, reads_file_for_segment_search.file, hits_for_read, juncs, read);
3301 }
3302 return result;
3303 }
3304
3305 static const int UNCOVERED = 0;
3306 static const int LOOK_LEFT = 1;
3307 static const int LOOK_RIGHT = 2;
3308
3309 uint8_t get_cov(const vector<uint8_t>& cov, uint32_t c)
3310 {
3311 uint32_t b = c >> 2;
3312 uint32_t r = c & 0x3;
3313 uint8_t s = (r << 1);
3314 uint8_t v = cov[b];
3315 v &= (0x3 << s);
3316 v >>= s;
3317 return v;
3318 }
3319
3320 void pair_covered_sites(ReadTable& it,
3321 RefSequenceTable& rt,
3322 vector<FZPipe*>& seg_files,
3323 std::set<Junction, skip_count_lt>& cov_juncs,
3324 size_t half_splice_mer_len)
3325 {
3326 BowtieHitFactory hit_factory(it,rt);
3327 map<uint32_t, vector<bool> > coverage_map;
3328 vector<RefSeg> expected_look_left_windows;
3329 vector<RefSeg> expected_look_right_windows;
3330
3331 for (size_t f = 0; f < seg_files.size(); ++f)
3332 {
3333 fprintf(stderr, "Adding hits from segment file %d to coverage map\n", (int)f);
3334 seg_files[f]->rewind();
3335 FILE* fp = seg_files[f]->file;
3336 //rewind(fp);
3337 HitStream hs(fp, &hit_factory, false, false, false);
3338
3339 HitsForRead hit_group;
3340 while (hs.next_read_hits(hit_group))
3341 {
3342 for (size_t h = 0; h < hit_group.hits.size(); ++h)
3343 {
3344 BowtieHit& bh = hit_group.hits[h];
3345
3346 pair<map<uint32_t, vector<bool> >::iterator, bool> ret =
3347 coverage_map.insert(make_pair(bh.ref_id(), vector<bool>()));
3348 vector<bool>& ref_cov = ret.first->second;
3349
3350 size_t right_extent = bh.right();
3351
3352 if (right_extent >= ref_cov.size())
3353 {
3354 ref_cov.resize(right_extent + 1, 0);
3355 }
3356
3357 for (uint32_t c = (uint32_t)bh.left(); c < (uint32_t)bh.right(); ++c)
3358 {
3359 ref_cov[c] = true;
3360 }
3361 }
3362 }
3363 }
3364
3365 static const int extend = 45;
3366 int num_islands = 0;
3367
3368 vector<RefSeg> expected_don_acc_windows;
3369
3370 fprintf(stderr, "Recording coverage islands\n");
3371 size_t cov_bases = 0;
3372 for (map<uint32_t, vector<bool> >::iterator itr = coverage_map.begin();
3373 itr != coverage_map.end();
3374 ++itr)
3375 {
3376 vector<bool>& cov = itr->second;
3377
3378 size_t island_left_edge = 0;
3379 for (size_t c = 1; c < cov.size(); ++c)
3380 {
3381 if (cov[c])
3382 {
3383 cov_bases++;
3384 if (!cov[c - 1])
3385 {
3386 num_islands += 1;
3387
3388 int edge = (int)c - extend;
3389 edge = max(edge, 0);
3390 island_left_edge = edge;
3391 }
3392 }
3393 else
3394 {
3395 if (cov[c - 1])
3396 {
3397 expected_don_acc_windows.push_back(RefSeg(itr->first,
3398 POINT_DIR_LEFT,
3399 false, /* not important */
3400 READ_DONTCARE,
3401 island_left_edge,
3402 c + extend));
3403 expected_don_acc_windows.push_back(RefSeg(itr->first,
3404 POINT_DIR_RIGHT,
3405 false, /* not important */
3406 READ_DONTCARE,
3407 island_left_edge,
3408 c + extend));
3409 }
3410 }
3411 }
3412 }
3413 fprintf(stderr, "Found %d islands covering %ld bases\n", num_islands, (long int)cov_bases);
3414
3415 juncs_from_ref_segs<RecordButterflyJuncs>(rt,
3416 expected_don_acc_windows,
3417 cov_juncs,
3418 "GT",
3419 "AG",
3420 max_coverage_intron_length,
3421 min_coverage_intron_length,
3422 max_cov_juncs,
3423 true,
3424 half_splice_mer_len);
3425 fprintf(stderr, "Found %ld potential intra-island junctions\n", (long int)cov_juncs.size());
3426 }
3427
3428 void capture_island_ends(ReadTable& it,
3429 RefSequenceTable& rt,
3430 //vector<FILE*>& seg_files,
3431 vector<FZPipe*>& seg_files,
3432 std::set<Junction, skip_count_lt>& cov_juncs,
3433 size_t half_splice_mer_len)
3434 {
3435 //static int island_repeat_tolerance = 10;
3436 BowtieHitFactory hit_factory(it,rt);
3437 map<uint32_t, vector<bool> > coverage_map;
3438 vector<RefSeg> expected_look_left_windows;
3439 vector<RefSeg> expected_look_right_windows;
3440
3441 for (size_t f = 0; f < seg_files.size(); ++f)
3442 {
3443 fprintf(stderr, "Adding hits from segment file %d to coverage map\n", (int)f);
3444 seg_files[f]->rewind();
3445 FILE* fp = seg_files[f]->file;
3446 //rewind(fp);
3447 HitStream hs(fp, &hit_factory, false, false, false);
3448
3449 HitsForRead hit_group;
3450 while (hs.next_read_hits(hit_group))
3451 {
3452 for (size_t h = 0; h < hit_group.hits.size(); ++h)
3453 {
3454 BowtieHit& bh = hit_group.hits[h];
3455
3456 pair<map<uint32_t, vector<bool> >::iterator, bool> ret =
3457 coverage_map.insert(make_pair(bh.ref_id(), vector<bool>()));
3458 vector<bool>& ref_cov = ret.first->second;
3459
3460 size_t right_extent = bh.right();
3461
3462 if (right_extent >= ref_cov.size())
3463 {
3464 ref_cov.resize(right_extent + 1, 0);
3465 }
3466
3467 for (uint32_t c = (uint32_t)bh.left(); c < (uint32_t)bh.right(); ++c)
3468 {
3469 ref_cov[c] = true;
3470 }
3471 }
3472 }
3473 }
3474
3475 static const int min_cov_length = segment_length + 2;
3476
3477 long covered_bases = 0;
3478 int long_enough_bases = 0;
3479 int left_looking = 0;
3480 int right_looking = 0;
3481 static const int extend = 45;
3482 static const int repeat_tol = 5;
3483
3484 int num_islands = 0;
3485
3486 for (map<uint32_t, vector<bool> >::iterator itr = coverage_map.begin();
3487 itr != coverage_map.end();
3488 ++itr)
3489 {
3490 //fprintf (stderr, "Finding pairings in ref seg %u\n", (int)itr->first);
3491
3492 vector<bool>& cov = itr->second;
3493 vector<bool> long_enough(cov.size());
3494
3495 size_t last_uncovered = 0;
3496
3497 static const uint8_t min_cov = 1;
3498
3499 for (size_t c = 1; c < cov.size(); ++c)
3500 {
3501 int c_cov = cov[c]; //get_cov(cov, c);
3502 if (c_cov < min_cov || c == (cov.size()) - 1)
3503 {
3504 int putative_exon_length = (int)c - (int)last_uncovered;
3505 int last_pos_cov = cov[c - 1]; //get_cov(cov,c - 1);
3506 if (last_pos_cov >= min_cov && putative_exon_length >= min_cov_length)
3507 {
3508 //fprintf (stdout, "%s\t%d\t%d\n", rt.get_name(itr->first), last_uncovered + 1, c);
3509 covered_bases += (c + 1 - last_uncovered);
3510 //fprintf(stderr, "making new segment from %d to %d\n", last_uncovered + 1, c);
3511 for (int l = (int)c; l > (int)last_uncovered; --l)
3512 {
3513 long_enough[l] = true;
3514 }
3515 }
3516 last_uncovered = c;
3517 }
3518 }
3519
3520 vector<bool>& ref_cov = long_enough;
3521 vector<uint8_t> cov_state(ref_cov.size(), UNCOVERED);
3522
3523 for (size_t c = 1; c < ref_cov.size(); ++c)
3524 {
3525 if (ref_cov[c])
3526 {
3527 long_enough_bases++;
3528 if (!ref_cov[c - 1])
3529 {
3530 num_islands += 1;
3531
3532 for (int r = (int)c - extend;
3533 r >= 0 && r < (int)c + repeat_tol && r < (int)cov_state.size();
3534 ++r)
3535 {
3536 cov_state[r] |= LOOK_LEFT;
3537 }
3538 }
3539 }
3540 else
3541 {
3542 if (ref_cov[c - 1])
3543 {
3544 for (int l = (int)c - repeat_tol;
3545 l >= 0 && l < (int)c + extend && l < (int)cov_state.size();
3546 ++l)
3547 {
3548 cov_state[l] |= LOOK_RIGHT;
3549 }
3550 }
3551 }
3552 }
3553
3554 RefSeg* curr_look_left = NULL;
3555 RefSeg* curr_look_right = NULL;
3556
3557 for (size_t c = 1; c < cov_state.size(); ++c)
3558 {
3559 if (cov_state[c] & LOOK_LEFT)
3560 {
3561 left_looking++;
3562 if (!(cov_state[c-1] & LOOK_LEFT))
3563 {
3564 expected_look_left_windows.push_back(RefSeg(itr->first,
3565 POINT_DIR_LEFT,
3566 false, /* not important */
3567 READ_DONTCARE,
3568 c,
3569 c + 1));
3570 curr_look_left = &(expected_look_left_windows.back());
3571 }
3572 else if (curr_look_left)
3573 {
3574 curr_look_left->right++;
3575 }
3576 }
3577 else
3578 {
3579 if ((cov_state[c-1] & LOOK_LEFT))
3580 {
3581 curr_look_left = NULL;
3582 }
3583 }
3584
3585 if (cov_state[c] & LOOK_RIGHT)
3586 {
3587 right_looking++;
3588 if (!(cov_state[c-1] & LOOK_RIGHT))
3589 {
3590 expected_look_right_windows.push_back(RefSeg(itr->first,
3591 POINT_DIR_RIGHT,
3592 false, /* not important */
3593 READ_DONTCARE,
3594 c,
3595 c + 1));
3596
3597 curr_look_right = &(expected_look_right_windows.back());
3598 }
3599 else if (curr_look_right)
3600 {
3601 curr_look_right->right++;
3602 }
3603 }
3604 else
3605 {
3606 if ((cov_state[c-1] & LOOK_RIGHT))
3607 {
3608 curr_look_right = NULL;
3609 }
3610 }
3611 }
3612 }
3613
3614 fprintf(stderr, "Map covers %ld bases\n", covered_bases);
3615 fprintf(stderr, "Map covers %d bases in sufficiently long segments\n", long_enough_bases);
3616 fprintf(stderr, "Map contains %d good islands\n", num_islands + 1);
3617 fprintf(stderr, "%d are left looking bases\n", left_looking);
3618 fprintf(stderr, "%d are right looking bases\n", right_looking);
3619
3620 vector<RefSeg> expected_don_acc_windows;
3621 expected_don_acc_windows.insert(expected_don_acc_windows.end(),
3622 expected_look_right_windows.begin(),
3623 expected_look_right_windows.end());
3624
3625 expected_don_acc_windows.insert(expected_don_acc_windows.end(),
3626 expected_look_left_windows.begin(),
3627 expected_look_left_windows.end());
3628
3629 coverage_map.clear();
3630 juncs_from_ref_segs<RecordExtendableJuncs>(rt,
3631 expected_don_acc_windows,
3632 cov_juncs,
3633 "GT",
3634 "AG",
3635 max_coverage_intron_length,
3636 min_coverage_intron_length,
3637 max_cov_juncs,
3638 true,
3639 half_splice_mer_len);
3640 fprintf(stderr, "Found %ld potential island-end pairing junctions\n", (long int)cov_juncs.size());
3641 }
3642
3643
3644 void process_segment_hits(RefSequenceTable& rt,
3645 //FILE* reads_file,
3646 FZPipe& reads_file,
3647 FZPipe& reads_file_for_segment_search,
3648 FZPipe& reads_file_for_indel_discovery,
3649 //vector<FILE*>& seg_files,
3650 vector<FZPipe>& seg_files,
3651 BowtieHitFactory& hit_factory,
3652 ReadTable& it,
3653 std::set<Junction, skip_count_lt>& juncs,
3654 std::set<Deletion>& deletions,
3655 std::set<Insertion>& insertions,
3656 eREAD read = READ_DONTCARE)
3657 {
3658 //get_seqs(reads_stream, it, false, false);
3659 //get_seqs(reads_file, it, reads_format, false, false);
3660
3661 vector<HitStream> hit_streams;
3662 for (size_t i = 0; i < seg_files.size(); ++i)
3663 {
3664 HitStream hs(seg_files[i].file, &hit_factory, false, false, false);
3665
3666 // if the next group id in this stream is zero, it's an empty stream,
3667 // and we can simply skip it.
3668 //if (hs.next_group_id() != 0)
3669 hit_streams.push_back(hs);
3670 }
3671
3672 int num_group = 0;
3673 while (process_next_hit_group(rt,
3674 reads_file,
3675 reads_file_for_segment_search,
3676 reads_file_for_indel_discovery,
3677 it,
3678 hit_streams,
3679 hit_streams.size() - 1,
3680 juncs, deletions, insertions, read) == true)
3681 {
3682 num_group++;
3683 if (num_group % 500000 == 0)
3684 fprintf(stderr, "\tProcessed %d root segment groups\n", num_group);
3685 }
3686 fprintf(stderr, "Microaligned %d segments\n", microaligned_segs);
3687 }
3688
3689 void print_juncs(RefSequenceTable& rt, std::set<Junction, skip_count_lt>& juncs, const char* str)
3690 {
3691 fprintf (stderr, "-- %s --\n", str);
3692 for(std::set<Junction>::iterator itr = juncs.begin();
3693 itr != juncs.end();
3694 ++itr)
3695 {
3696 const char* ref_name = rt.get_name(itr->refid);
3697
3698 fprintf(stderr, "%s\t%d\t%d\t%c\n",
3699 ref_name,
3700 itr->left,
3701 itr->right,
3702 itr->antisense ? '-' : '+');
3703 }
3704 fprintf (stderr, "-- done --\n");
3705 }
3706
3707 void driver(istream& ref_stream,
3708 FILE* juncs_out,
3709 FILE* insertions_out,
3710 FILE* deletions_out,
3711 //FILE* left_reads_file,
3712 FZPipe& left_reads_file,
3713 FZPipe& left_reads_file_for_segment_search,
3714 FZPipe& left_reads_file_for_indel_discovery,
3715 vector<FZPipe>& left_seg_files,
3716 FZPipe& right_reads_file,
3717 FZPipe& right_reads_file_for_segment_search,
3718 FZPipe& right_reads_file_for_indel_discovery,
3719 vector<FZPipe>& right_seg_files)
3720 {
3721 if (left_seg_files.size() == 0)
3722 {
3723 fprintf(stderr, "No hits to process, exiting\n");
3724 exit(0);
3725 }
3726
3727 std::set<Junction, skip_count_lt> seg_juncs;
3728 std::set<Junction, skip_count_lt> cov_juncs;
3729 std::set<Junction, skip_count_lt> butterfly_juncs;
3730
3731 std::set<Junction> juncs;
3732
3733 std::set<Deletion> deletions;
3734 std::set<Insertion> insertions;
3735
3736
3737 RefSequenceTable rt(true, true);
3738
3739 fprintf (stderr, "Loading reference sequences...\n");
3740 get_seqs(ref_stream, rt, true, false);
3741
3742 ReadTable it;
3743 BowtieHitFactory hit_factory(it,rt);
3744
3745 if (left_seg_files.size() > 1)
3746 {
3747 fprintf( stderr, "Loading left segment hits...\n");
3748 process_segment_hits(rt,
3749 left_reads_file,
3750 left_reads_file_for_segment_search,
3751 left_reads_file_for_indel_discovery,
3752 left_seg_files,
3753 hit_factory,
3754 it,
3755 seg_juncs,
3756 deletions,
3757 insertions,
3758 READ_LEFT);
3759 fprintf( stderr, "done.\n");
3760 }
3761
3762 if (right_seg_files.size() > 1)
3763 {
3764 fprintf( stderr, "Loading right segment hits...\n");
3765 process_segment_hits(rt,
3766 right_reads_file,
3767 right_reads_file_for_segment_search,
3768 right_reads_file_for_indel_discovery,
3769 right_seg_files,
3770 hit_factory,
3771 it,
3772 seg_juncs,
3773 deletions,
3774 insertions,
3775 READ_RIGHT);
3776 fprintf( stderr, "done.\n");
3777 }
3778 fprintf(stderr, "Found %ld potential split-segment junctions\n", (long int)seg_juncs.size());
3779 fprintf(stderr, "Found %ld potential small deletions\n", (long int)deletions.size());
3780 fprintf(stderr, "Found %ld potential small insertions\n", (long int)insertions.size());
3781
3782 //vector<FILE*> all_seg_files;
3783 vector<FZPipe*> all_seg_files;
3784 for (vector<FZPipe>::size_type i = 0; i != left_seg_files.size();i++) {
3785 all_seg_files.push_back( &(left_seg_files[i]) );
3786 }
3787 for (vector<FZPipe>::size_type i = 0; i != right_seg_files.size();i++) {
3788 all_seg_files.push_back( &(right_seg_files[i]) );
3789 }
3790 //copy(left_seg_files.begin(), left_seg_files.end(), back_inserter(all_seg_files));
3791 //copy(right_seg_files.begin(), right_seg_files.end(), back_inserter(all_seg_files));
3792
3793 #if 0
3794 // daehwan - check this out as Cole insists on using segments gives better results.
3795 vector<FILE*> all_map_files;
3796 all_map_files.push_back(left_reads_map_file);
3797 all_map_files.push_back(right_reads_map_file);
3798 copy(all_seg_files.begin(), all_seg_files.end(), back_inserter(all_map_files));
3799 #endif
3800
3801 if (!no_coverage_search || butterfly_search)
3802 {
3803 if (ium_reads != "")
3804 {
3805 vector<string> ium_read_files;
3806 tokenize(ium_reads,",", ium_read_files);
3807 //vector<FILE*> iums;
3808 vector<FZPipe> iums;
3809 string unzcmd=getUnpackCmd(ium_read_files[0],false);
3810 for (size_t ium = 0; ium < ium_read_files.size(); ++ium)
3811 {
3812 fprintf (stderr, "Indexing extensions in %s\n", ium_read_files[ium].c_str());
3813 //FILE* ium_file = fopen(ium_read_files[ium].c_str(), "r");
3814 FZPipe ium_file(ium_read_files[ium],unzcmd);
3815 if (ium_file.file==NULL)
3816 {
3817 fprintf (stderr, "Can't open file %s for reading, skipping...\n",ium_read_files[ium].c_str());
3818 continue;
3819 }
3820 iums.push_back(ium_file);
3821 }
3822
3823 index_read_mers(iums, 5);
3824 }
3825 else
3826 {
3827 no_coverage_search = true;
3828 butterfly_search = false;
3829 }
3830
3831 if (!no_coverage_search)
3832 {
3833 fprintf(stderr, "Looking for junctions by island end pairings\n");
3834 capture_island_ends(it,
3835 rt,
3836 all_seg_files,
3837 cov_juncs,
3838 5);
3839 fprintf(stderr, "done\n");
3840 }
3841 }
3842
3843 if (butterfly_search)
3844 {
3845 fprintf(stderr, "Looking for junctions between and within islands\n");
3846 prune_extension_table(butterfly_overhang);
3847 compact_extension_table();
3848 pair_covered_sites(it,
3849 rt,
3850 all_seg_files,
3851 butterfly_juncs,
3852 5);
3853
3854 fprintf(stderr, "done\n");
3855 }
3856
3857 std::set<Junction, skip_count_lt> microexon_juncs;
3858 if (!no_microexon_search)
3859 {
3860 std::set<Junction, skip_count_lt> microexon_juncs;
3861 align_microexon_segs(rt,
3862 microexon_juncs,
3863 max_cov_juncs,
3864 5);
3865 juncs.insert(microexon_juncs.begin(), microexon_juncs.end());
3866 }
3867
3868 juncs.insert(cov_juncs.begin(), cov_juncs.end());
3869 juncs.insert(seg_juncs.begin(), seg_juncs.end());
3870 juncs.insert(butterfly_juncs.begin(), butterfly_juncs.end());
3871
3872 fprintf(stderr, "Reporting potential splice junctions...");
3873 for(std::set<Junction>::iterator itr = juncs.begin();
3874 itr != juncs.end();
3875 ++itr)
3876 {
3877 const char* ref_name = rt.get_name(itr->refid);
3878
3879 fprintf(juncs_out,
3880 "%s\t%d\t%d\t%c\n",
3881 ref_name,
3882 itr->left,
3883 itr->right,
3884 itr->antisense ? '-' : '+');
3885 }
3886 //close all reading pipes, just to exit cleanly
3887
3888 for (vector<FZPipe*>::size_type i = 0; i != all_seg_files.size();i++) {
3889 all_seg_files[i]->close();
3890 }
3891 left_reads_file.close();
3892 left_reads_file_for_indel_discovery.close();
3893 right_reads_file.close();
3894 right_reads_file_for_indel_discovery.close();
3895 fprintf (stderr, "done\n");
3896 fprintf (stderr, "Reported %d total possible splices\n", (int)juncs.size());
3897
3898
3899 fprintf (stderr, "Reporting potential deletions...\n");
3900 if(deletions_out){
3901 for(std::set<Deletion>::iterator itr = deletions.begin(); itr != deletions.end(); ++itr){
3902 const char* ref_name = rt.get_name(itr->refid);
3903 /*
3904 * We fix up the left co-ordinate to reference the first deleted base
3905 */
3906 fprintf(deletions_out,
3907 "%s\t%d\t%d\n",
3908 ref_name,
3909 itr->left + 1,
3910 itr->right);
3911 }
3912 fclose(deletions_out);
3913 }else{
3914 fprintf(stderr, "Failed to open deletions file for writing, no deletions reported\n");
3915 }
3916
3917 fprintf (stderr, "Reporting potential insertions...\n");
3918 if(insertions_out){
3919 for(std::set<Insertion>::iterator itr = insertions.begin(); itr != insertions.end(); ++itr){
3920 const char* ref_name = rt.get_name(itr->refid);
3921 fprintf(insertions_out,
3922 "%s\t%d\t%d\t%s\n",
3923 ref_name,
3924 itr->left,
3925 itr->left,
3926 itr->sequence.c_str());
3927 }
3928 fclose(insertions_out);
3929 }else{
3930 fprintf(stderr, "Failed to open insertions file for writing, no insertions reported\n");
3931 }
3932 }
3933
3934 int main(int argc, char** argv)
3935 {
3936 fprintf(stderr, "segment_juncs v%s (%s)\n", PACKAGE_VERSION, SVN_REVISION);
3937 fprintf(stderr, "---------------------------\n");
3938
3939 int parse_ret = parse_options(argc, argv, print_usage);
3940 if (parse_ret)
3941 return parse_ret;
3942
3943 if(optind >= argc)
3944 {
3945 print_usage();
3946 return 1;
3947 }
3948
3949 string ref_file_name = argv[optind++];
3950
3951 if(optind >= argc)
3952 {
3953 print_usage();
3954 return 1;
3955 }
3956
3957 string juncs_file_name = argv[optind++];
3958
3959 if(optind >= argc)
3960 {
3961 print_usage();
3962 return 1;
3963 }
3964
3965 string insertions_file_name = argv[optind++];
3966 if(optind >= argc)
3967 {
3968 print_usage();
3969 return 1;
3970 }
3971
3972 string deletions_file_name = argv[optind++];
3973 if(optind >= argc)
3974 {
3975 print_usage();
3976 return 1;
3977 }
3978
3979 string left_reads_file_name = argv[optind++];
3980
3981 if(optind >= argc)
3982 {
3983 print_usage();
3984 return 1;
3985 }
3986
3987 string left_reads_map_file_name = argv[optind++];
3988
3989 if(optind >= argc)
3990 {
3991 print_usage();
3992 return 1;
3993 }
3994
3995 string left_segment_file_list = argv[optind++];
3996
3997 string right_segment_file_list;
3998 string right_reads_file_name;
3999 string right_reads_map_file_name;
4000 if (optind < argc)
4001 {
4002 right_reads_file_name = argv[optind++];
4003
4004 if(optind >= argc)
4005 {
4006 print_usage();
4007 return 1;
4008 }
4009
4010 right_reads_map_file_name = argv[optind++];
4011
4012 if(optind >= argc)
4013 {
4014 print_usage();
4015 return 1;
4016 }
4017 right_segment_file_list = argv[optind++];
4018 }
4019
4020 // Open the approppriate files
4021
4022 ifstream ref_stream(ref_file_name.c_str(), ifstream::in);
4023 if (!ref_stream.good())
4024 {
4025 fprintf(stderr, "Error: cannot open %s for reading\n",
4026 ref_file_name.c_str());
4027 exit(1);
4028 }
4029
4030
4031 FILE* juncs_file = fopen(juncs_file_name.c_str(), "w");
4032 if (!juncs_file)
4033 {
4034 fprintf(stderr, "Error: cannot open %s for writing\n",
4035 juncs_file_name.c_str());
4036 exit(1);
4037 }
4038
4039 FILE* insertions_file = fopen(insertions_file_name.c_str(), "w");
4040 if (!insertions_file)
4041 {
4042 fprintf(stderr, "Error: cannot open %s for writing\n",
4043 insertions_file_name.c_str());
4044 exit(1);
4045 }
4046
4047
4048 FILE* deletions_file = fopen(deletions_file_name.c_str(), "w");
4049 if (!deletions_file)
4050 {
4051 fprintf(stderr, "Error: cannot open %s for writing\n",
4052 deletions_file_name.c_str());
4053 exit(1);
4054 }
4055
4056 //FILE* left_reads_file = fopen(left_reads_file_name.c_str(), "r");
4057 //FILE* left_reads_file_for_indel_discovery = fopen(left_reads_file_name.c_str(),"r");
4058 string unzcmd=getUnpackCmd(left_reads_file_name, false);
4059 FZPipe left_reads_file(left_reads_file_name, unzcmd);
4060 FZPipe left_reads_file_for_segment_search(left_reads_file_name, unzcmd);
4061 FZPipe left_reads_file_for_indel_discovery(left_reads_file_name, unzcmd);
4062 if (left_reads_file.file==NULL || left_reads_file_for_segment_search.file==NULL ||
4063 left_reads_file_for_indel_discovery.file==NULL)
4064 {
4065 fprintf(stderr, "Error: cannot open %s for reading\n",
4066 left_reads_file_name.c_str());
4067 exit(1);
4068 }
4069
4070 //FILE* left_reads_map_file = fopen(left_reads_map_file_name.c_str(), "r");
4071 FZPipe left_reads_map_file(left_reads_map_file_name, unzcmd);
4072 if (left_reads_map_file.file==NULL)
4073 {
4074 fprintf(stderr, "Error: cannot open %s for reading\n",
4075 left_reads_map_file_name.c_str());
4076 exit(1);
4077 }
4078
4079 vector<string> left_segment_file_names;
4080 //vector<FILE*> left_segment_files;
4081 vector<FZPipe> left_segment_files;
4082 tokenize(left_segment_file_list, ",",left_segment_file_names);
4083 for (size_t i = 0; i < left_segment_file_names.size(); ++i)
4084 {
4085 //FILE* seg_file = fopen(left_segment_file_names[i].c_str(), "r");
4086 FZPipe seg_file(left_segment_file_names[i], unzcmd);
4087 if (seg_file.file == NULL)
4088 {
4089 fprintf(stderr, "Error: cannot open segment map file %s for reading\n",
4090 left_segment_file_names[i].c_str());
4091 exit(1);
4092 }
4093 left_segment_files.push_back(seg_file);
4094 }
4095
4096 //vector<FILE*> right_segment_files;
4097 vector<FZPipe> right_segment_files;
4098 //FILE* right_reads_file = NULL;
4099 FZPipe right_reads_file;
4100 //FILE* right_reads_file_for_indel_discovery = NULL;
4101 FZPipe right_reads_file_for_segment_search;
4102 FZPipe right_reads_file_for_indel_discovery;
4103
4104 if (right_segment_file_list != "")
4105 {
4106 //right_reads_file = fopen(right_reads_file_name.c_str(), "r");
4107 right_reads_file.openRead(right_reads_file_name, unzcmd);
4108 //right_reads_file_for_indel_discovery = fopen(right_reads_file_name.c_str(), "r");
4109 right_reads_file_for_segment_search.openRead(right_reads_file_name, unzcmd);
4110 right_reads_file_for_indel_discovery.openRead(right_reads_file_name, unzcmd);
4111 if (right_reads_file.file==NULL || right_reads_file_for_indel_discovery.file==NULL)
4112 {
4113 fprintf(stderr, "Error: cannot open %s for reading\n",
4114 right_reads_file_name.c_str());
4115 exit(1);
4116 }
4117
4118 vector<string> right_segment_file_names;
4119
4120 tokenize(right_segment_file_list, ",",right_segment_file_names);
4121 for (size_t i = 0; i < right_segment_file_names.size(); ++i)
4122 {
4123 //FILE* seg_file = fopen(right_segment_file_names[i].c_str(), "r");
4124 FZPipe seg_file(right_segment_file_names[i],unzcmd);
4125 if (seg_file.file == NULL)
4126 {
4127 fprintf(stderr, "Error: cannot open segment map %s for reading\n",
4128 right_segment_file_names[i].c_str());
4129 exit(1);
4130 }
4131 right_segment_files.push_back(seg_file);
4132 }
4133 }
4134
4135 driver(ref_stream,
4136 juncs_file,
4137 insertions_file,
4138 deletions_file,
4139 left_reads_file,
4140 left_reads_file_for_segment_search,
4141 left_reads_file_for_indel_discovery,
4142 left_segment_files,
4143 right_reads_file,
4144 right_reads_file_for_segment_search,
4145 right_reads_file_for_indel_discovery,
4146 right_segment_files);
4147
4148 return 0;
4149 }