31 |
|
#include <getopt.h> |
32 |
|
|
33 |
|
#include <boost/thread.hpp> |
34 |
+ |
#include <boost/random/mersenne_twister.hpp> |
35 |
|
|
36 |
|
#include "common.h" |
37 |
|
#include "utils.h" |
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 |
|
|
82 |
|
rt.get_id(name, keep_seqs ? ref_str : NULL, 0); |
83 |
|
} |
147 |
|
} |
148 |
|
} |
149 |
|
|
150 |
< |
bool set_insert_alignment_grade(const BowtieHit& lh, const BowtieHit& rh, InsertAlignmentGrade& grade) |
150 |
> |
bool set_insert_alignment_grade(const BowtieHit& lh, const BowtieHit& rh, const JunctionSet& junctions, InsertAlignmentGrade& grade) |
151 |
|
{ |
149 |
– |
// max mate inner distance (genomic) |
150 |
– |
int min_mate_inner_dist = inner_dist_mean - inner_dist_std_dev; |
151 |
– |
if (max_mate_inner_dist == -1) |
152 |
– |
max_mate_inner_dist = inner_dist_mean + inner_dist_std_dev; |
153 |
– |
|
152 |
|
bool fusion = false; |
153 |
|
bool left_fusion = lh.fusion_opcode() != FUSION_NOTHING; |
154 |
|
bool right_fusion = rh.fusion_opcode() != FUSION_NOTHING; |
178 |
|
} |
179 |
|
} |
180 |
|
|
181 |
< |
grade = InsertAlignmentGrade(lh, rh, min_mate_inner_dist, max_mate_inner_dist, fusion); |
181 |
> |
grade = InsertAlignmentGrade(lh, rh, junctions, fusion); |
182 |
|
|
183 |
|
return true; |
184 |
|
} |
185 |
|
|
186 |
|
struct cmp_pair_alignment |
187 |
|
{ |
188 |
< |
cmp_pair_alignment(const JunctionSet& gtf_juncs) : |
189 |
< |
_gtf_juncs(>f_juncs) {} |
188 |
> |
cmp_pair_alignment(const JunctionSet& junctions) : |
189 |
> |
_junctions(&junctions) {} |
190 |
|
|
191 |
|
bool operator() (const pair<BowtieHit, BowtieHit>& l, const pair<BowtieHit, BowtieHit>& r) const |
192 |
|
{ |
193 |
< |
InsertAlignmentGrade gl; set_insert_alignment_grade(l.first, l.second, gl); |
194 |
< |
InsertAlignmentGrade gr; set_insert_alignment_grade(r.first, r.second, gr); |
193 |
> |
InsertAlignmentGrade gl; set_insert_alignment_grade(l.first, l.second, *_junctions, gl); |
194 |
> |
InsertAlignmentGrade gr; set_insert_alignment_grade(r.first, r.second, *_junctions, gr); |
195 |
|
|
196 |
|
bool better = gr < gl; |
197 |
|
bool worse = gl < gr; |
202 |
|
return false; |
203 |
|
} |
204 |
|
|
205 |
< |
const JunctionSet* _gtf_juncs; |
205 |
> |
const JunctionSet* _junctions; |
206 |
|
}; |
207 |
|
|
208 |
|
void pair_best_alignments(const HitsForRead& left_hits, |
241 |
|
junctions, insertions, deletions, fusions, coverage); |
242 |
|
rh.alignment_score(align_status._alignment_score); |
243 |
|
InsertAlignmentGrade g; |
244 |
< |
bool allowed = set_insert_alignment_grade(lh, rh, g); |
244 |
> |
bool allowed; |
245 |
> |
allowed = set_insert_alignment_grade(lh, rh, final_report ? junctions : gtf_junctions, g); |
246 |
> |
|
247 |
> |
// daehwan - for debugging purposes |
248 |
> |
if (lh.insert_id() == 325708 && !g.is_fusion() && false) |
249 |
> |
{ |
250 |
> |
fprintf(stderr, "lh %d:%d %s score: %d (from %d) NM: %d\n", |
251 |
> |
lh.ref_id(), lh.left(), print_cigar(lh.cigar()).c_str(), |
252 |
> |
lh.alignment_score(), left[i].alignment_score(), lh.edit_dist()); |
253 |
> |
fprintf(stderr, "rh %d:%d %s score: %d (from %d) NM: %d\n", |
254 |
> |
rh.ref_id(), rh.left(), print_cigar(rh.cigar()).c_str(), |
255 |
> |
rh.alignment_score(), right[j].alignment_score(), rh.edit_dist()); |
256 |
> |
fprintf(stderr, "combined score: %d is_fusion(%d)\n", g.align_score(), g.is_fusion()); |
257 |
> |
} |
258 |
|
|
259 |
|
if (!allowed) continue; |
260 |
|
if (!fusion_search && !report_discordant_pair_alignments && g.is_fusion()) continue; |
282 |
|
|
283 |
|
if ((report_secondary_alignments || !final_report) && best_hits.size() > 0) |
284 |
|
{ |
285 |
< |
cmp_pair_alignment cmp(gtf_junctions); |
285 |
> |
cmp_pair_alignment cmp(final_report ? junctions : gtf_junctions); |
286 |
|
sort(best_hits.begin(), best_hits.end(), cmp); |
287 |
< |
set_insert_alignment_grade(best_hits[0].first, best_hits[0].second, best_grade); |
287 |
> |
set_insert_alignment_grade(best_hits[0].first, best_hits[0].second, final_report ? junctions : gtf_junctions, best_grade); |
288 |
|
} |
289 |
|
|
290 |
|
if (final_report) |
393 |
|
string ref_name = sam_toks[2], ref_name2 = ""; |
394 |
|
char cigar1[1024] = {0}, cigar2[1024] = {0}; |
395 |
|
string left_seq, right_seq, left_qual, right_qual; |
396 |
< |
int left = -1, left1 = -1, left2 = -1; |
396 |
> |
int left1 = -1, left2 = -1; |
397 |
|
bool fusion_alignment = false; |
398 |
|
size_t XF_index = 0; |
399 |
|
for (size_t i = 11; i < sam_toks.size(); ++i) |
424 |
|
else if (strncmp(tok.c_str(), "AS", 2) == 0) |
425 |
|
{ |
426 |
|
char AS_score[128] = {0}; |
427 |
< |
sprintf(AS_score, "AS:i:%d", bh.alignment_score()); |
427 |
> |
sprintf(AS_score, "AS:i:%d", min(0, bh.alignment_score())); |
428 |
|
tok = AS_score; |
429 |
|
} |
430 |
|
} |
452 |
|
flag |= 0x100; |
453 |
|
|
454 |
|
int gpos=isdigit(sam_toks[3][0]) ? atoi(sam_toks[3].c_str()) : 0; |
455 |
< |
int mapQ=255; |
455 |
> |
int mapQ = 50; |
456 |
|
if (num_hits > 1) { |
457 |
|
double err_prob = 1 - (1.0 / num_hits); |
458 |
|
mapQ = (int)(-10.0 * log(err_prob) / log(10.0)); |
460 |
|
int tlen =atoi(sam_toks[8].c_str()); //TLEN |
461 |
|
int mate_pos=atoi(sam_toks[7].c_str()); |
462 |
|
|
463 |
+ |
string rg_aux = ""; |
464 |
+ |
if (!sam_readgroup_id.empty()) |
465 |
+ |
rg_aux = string("RG:Z:") + sam_readgroup_id; |
466 |
+ |
|
467 |
|
GBamRecord* bamrec=NULL; |
468 |
|
if (fusion_alignment) { |
469 |
|
vector<string> auxdata; |
470 |
|
add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex); |
471 |
+ |
if (rg_aux != "") |
472 |
+ |
auxdata.push_back(rg_aux); |
473 |
|
bamrec=bam_writer.new_record(qname.c_str(), flag, ref_name.c_str(), left1 + 1, mapQ, |
474 |
|
cigar1, sam_toks[6].c_str(), mate_pos, |
475 |
|
tlen, left_seq.c_str(), left_qual.c_str(), &auxdata); |
479 |
|
auxdata.clear(); |
480 |
|
sam_toks[XF_index][5] = '2'; |
481 |
|
add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex); |
482 |
+ |
if (rg_aux != "") |
483 |
+ |
auxdata.push_back(rg_aux); |
484 |
|
bamrec=bam_writer.new_record(qname.c_str(), flag, ref_name2.c_str(), left2 + 1, mapQ, |
485 |
|
cigar2, sam_toks[6].c_str(), mate_pos, |
486 |
|
tlen, right_seq.c_str(), right_qual.c_str(), &auxdata); |
489 |
|
} else { |
490 |
|
vector<string> auxdata; |
491 |
|
add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex); |
492 |
+ |
if (rg_aux != "") |
493 |
+ |
auxdata.push_back(rg_aux); |
494 |
|
bamrec=bam_writer.new_record(qname.c_str(), flag, sam_toks[2].c_str(), gpos, mapQ, |
495 |
|
sam_toks[5].c_str(), sam_toks[6].c_str(), mate_pos, |
496 |
|
tlen, sam_toks[9].c_str(), sam_toks[10].c_str(), &auxdata); |
537 |
|
string ref_name = sam_toks[2], ref_name2 = ""; |
538 |
|
char cigar1[1024] = {0}, cigar2[1024] = {0}; |
539 |
|
string left_seq, right_seq, left_qual, right_qual; |
540 |
< |
int left = -1, left1 = -1, left2 = -1; |
540 |
> |
int left1 = -1, left2 = -1; |
541 |
|
bool fusion_alignment = false; |
542 |
|
size_t XF_tok_idx = 11; |
543 |
|
for (; XF_tok_idx < sam_toks.size(); ++XF_tok_idx) |
567 |
|
else if (strncmp(tok.c_str(), "AS", 2) == 0) |
568 |
|
{ |
569 |
|
char AS_score[128] = {0}; |
570 |
< |
sprintf(AS_score, "AS:i:%d", bh.alignment_score()); |
570 |
> |
sprintf(AS_score, "AS:i:%d", min(0, bh.alignment_score())); |
571 |
|
tok = AS_score; |
572 |
|
} |
573 |
|
} |
574 |
|
|
575 |
|
int gpos=isdigit(sam_toks[3][0]) ? atoi(sam_toks[3].c_str()) : 0; |
576 |
< |
int mapQ=255; |
576 |
> |
int mapQ = 50; |
577 |
|
if (grade.num_alignments > 1) { |
578 |
|
double err_prob = 1 - (1.0 / grade.num_alignments); |
579 |
|
mapQ = (int)(-10.0 * log(err_prob) / log(10.0)); |
610 |
|
flag |= 0x0008; |
611 |
|
} |
612 |
|
|
613 |
+ |
string rg_aux = ""; |
614 |
+ |
if (!sam_readgroup_id.empty()) |
615 |
+ |
rg_aux = string("RG:Z:") + sam_readgroup_id; |
616 |
+ |
|
617 |
|
GBamRecord* bamrec=NULL; |
618 |
|
if (fusion_alignment) { |
619 |
|
vector<string> auxdata; |
620 |
|
add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex); |
621 |
+ |
if (rg_aux != "") |
622 |
+ |
auxdata.push_back(rg_aux); |
623 |
|
bamrec=bam_writer.new_record(qname.c_str(), flag, ref_name.c_str(), left1 + 1, mapQ, |
624 |
|
cigar1, sam_toks[6].c_str(), mate_pos, |
625 |
|
tlen, left_seq.c_str(), left_qual.c_str(), &auxdata); |
629 |
|
auxdata.clear(); |
630 |
|
sam_toks[XF_tok_idx][5] = '2'; |
631 |
|
add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex); |
632 |
+ |
if (rg_aux != "") |
633 |
+ |
auxdata.push_back(rg_aux); |
634 |
|
bamrec=bam_writer.new_record(qname.c_str(), flag, ref_name2.c_str(), left2 + 1, mapQ, |
635 |
|
cigar2, sam_toks[6].c_str(), mate_pos, |
636 |
|
tlen, right_seq.c_str(), right_qual.c_str(), &auxdata); |
639 |
|
} else { |
640 |
|
vector<string> auxdata; |
641 |
|
add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex); |
642 |
+ |
if (rg_aux != "") |
643 |
+ |
auxdata.push_back(rg_aux); |
644 |
|
bamrec=bam_writer.new_record(qname.c_str(), flag, sam_toks[2].c_str(), gpos, mapQ, |
645 |
|
sam_toks[5].c_str(), sam_toks[6].c_str(), mate_pos, |
646 |
|
tlen, sam_toks[9].c_str(), sam_toks[10].c_str(), &auxdata); |
677 |
|
void print_sam_for_single(const RefSequenceTable& rt, |
678 |
|
const HitsForRead& hits, |
679 |
|
FragmentType frag_type, |
680 |
< |
Read& read, |
680 |
> |
const string& alt_name, |
681 |
|
GBamWriter& bam_writer, |
682 |
< |
GBamWriter* um_out //write unmapped reads here |
652 |
< |
) |
682 |
> |
boost::random::mt19937& rng) |
683 |
|
{ |
684 |
< |
assert(!read.alt_name.empty()); |
684 |
> |
//assert(!read.alt_name.empty()); |
685 |
|
if (hits.hits.empty()) |
686 |
|
return; |
687 |
|
|
696 |
|
|
697 |
|
size_t primaryHit = 0; |
698 |
|
if (!report_secondary_alignments) |
699 |
< |
primaryHit = random() % hits.hits.size(); |
699 |
> |
primaryHit = rng() % hits.hits.size(); |
700 |
|
|
701 |
|
bool multipleHits = (hits.hits.size() > 1); |
702 |
|
for (i = 0; i < hits.hits.size(); ++i) |
707 |
|
rewrite_sam_record(bam_writer, rt, |
708 |
|
bh, |
709 |
|
bh.hitfile_rec().c_str(), |
710 |
< |
read.alt_name.c_str(), |
710 |
> |
alt_name.c_str(), |
711 |
|
frag_type, |
712 |
|
hits.hits.size(), |
713 |
|
(i < hits.hits.size()-1) ? &(hits.hits[index_vector[i+1]]) : NULL, |
719 |
|
void print_sam_for_pair(const RefSequenceTable& rt, |
720 |
|
const vector<pair<BowtieHit, BowtieHit> >& best_hits, |
721 |
|
const InsertAlignmentGrade& grade, |
692 |
– |
ReadStream& left_reads_file, |
693 |
– |
ReadStream& right_reads_file, |
722 |
|
GBamWriter& bam_writer, |
723 |
< |
GBamWriter* left_um_out, |
724 |
< |
GBamWriter* right_um_out, |
723 |
> |
const string& left_alt_name, |
724 |
> |
const string& right_alt_name, |
725 |
> |
boost::random::mt19937& rng, |
726 |
|
uint64_t begin_id = 0, |
727 |
|
uint64_t end_id = std::numeric_limits<uint64_t>::max()) |
728 |
|
{ |
745 |
|
sort(index_vector.begin(), index_vector.end(), s); |
746 |
|
|
747 |
|
if (!report_secondary_alignments) |
748 |
< |
primaryHit = random() % right_hits.hits.size(); |
748 |
> |
primaryHit = rng() % right_hits.hits.size(); |
749 |
|
|
721 |
– |
bool got_left_read = left_reads_file.getRead(best_hits[0].first.insert_id(), left_read, |
722 |
– |
reads_format, false, begin_id, end_id, |
723 |
– |
left_um_out, false); |
724 |
– |
|
725 |
– |
bool got_right_read = right_reads_file.getRead(best_hits[0].first.insert_id(), right_read, |
726 |
– |
reads_format, false, begin_id, end_id, |
727 |
– |
right_um_out, false); |
728 |
– |
|
729 |
– |
assert (got_left_read && got_right_read); |
750 |
|
bool multipleHits = (best_hits.size() > 1); |
751 |
|
for (i = 0; i < best_hits.size(); ++i) |
752 |
|
{ |
758 |
|
rewrite_sam_record(bam_writer, rt, |
759 |
|
right_bh, |
760 |
|
right_bh.hitfile_rec().c_str(), |
761 |
< |
right_read.alt_name.c_str(), |
761 |
> |
right_alt_name.c_str(), |
762 |
|
grade, |
763 |
|
FRAG_RIGHT, |
764 |
|
&left_bh, |
769 |
|
rewrite_sam_record(bam_writer, rt, |
770 |
|
left_bh, |
771 |
|
left_bh.hitfile_rec().c_str(), |
772 |
< |
left_read.alt_name.c_str(), |
772 |
> |
left_alt_name.c_str(), |
773 |
|
grade, |
774 |
|
FRAG_LEFT, |
775 |
|
&right_bh, |
918 |
|
hits = remaining; |
919 |
|
} |
920 |
|
|
921 |
+ |
void realign_reads(HitsForRead& hits, |
922 |
+ |
const RefSequenceTable& rt, |
923 |
+ |
const JunctionSet& junctions, |
924 |
+ |
const JunctionSet& rev_junctions, |
925 |
+ |
const InsertionSet& insertions, |
926 |
+ |
const DeletionSet& deletions, |
927 |
+ |
const DeletionSet& rev_deletions, |
928 |
+ |
const FusionSet& fusions) |
929 |
+ |
{ |
930 |
+ |
vector<BowtieHit> additional_hits; |
931 |
+ |
for (size_t i = 0; i < hits.hits.size(); ++i) |
932 |
+ |
{ |
933 |
+ |
BowtieHit& bh = hits.hits[i]; |
934 |
+ |
|
935 |
+ |
const vector<CigarOp>& cigars = bh.cigar(); |
936 |
+ |
int pos = bh.left(); |
937 |
+ |
int refid = bh.ref_id(); |
938 |
+ |
for (size_t j = 0; j < cigars.size(); ++j) |
939 |
+ |
{ |
940 |
+ |
const CigarOp& op = cigars[j]; |
941 |
+ |
if (j == 0 || j == cigars.size() - 1) |
942 |
+ |
{ |
943 |
+ |
// let's do this for MATCH case only, |
944 |
+ |
if (op.opcode == MATCH || op.opcode == mATCH) |
945 |
+ |
{ |
946 |
+ |
int left1, left2; |
947 |
+ |
if (op.opcode == MATCH) |
948 |
+ |
{ |
949 |
+ |
left1 = pos; |
950 |
+ |
left2 = pos + op.length - 1; |
951 |
+ |
} |
952 |
+ |
else |
953 |
+ |
{ |
954 |
+ |
left1 = pos - op.length + 1; |
955 |
+ |
left2 = pos; |
956 |
+ |
} |
957 |
+ |
|
958 |
+ |
{ |
959 |
+ |
JunctionSet::const_iterator lb, ub; |
960 |
+ |
JunctionSet temp_junctions; |
961 |
+ |
if (j == 0) |
962 |
+ |
{ |
963 |
+ |
lb = rev_junctions.upper_bound(Junction(refid, left1, 0, true)); |
964 |
+ |
ub = rev_junctions.lower_bound(Junction(refid, left2, left2, false)); |
965 |
+ |
while (lb != ub && lb != rev_junctions.end()) |
966 |
+ |
{ |
967 |
+ |
Junction temp_junction = lb->first; |
968 |
+ |
temp_junction.left = lb->first.right; |
969 |
+ |
temp_junction.right = lb->first.left; |
970 |
+ |
temp_junctions[temp_junction] = lb->second; |
971 |
+ |
++lb; |
972 |
+ |
} |
973 |
+ |
} |
974 |
+ |
|
975 |
+ |
if (j == cigars.size() - 1) |
976 |
+ |
{ |
977 |
+ |
int common_right = left2 + max_report_intron_length; |
978 |
+ |
lb = junctions.upper_bound(Junction(refid, left1, common_right, true)); |
979 |
+ |
ub = junctions.lower_bound(Junction(refid, left2, common_right, false)); |
980 |
+ |
while (lb != ub && lb != junctions.end()) |
981 |
+ |
{ |
982 |
+ |
temp_junctions[lb->first] = lb->second; |
983 |
+ |
++lb; |
984 |
+ |
} |
985 |
+ |
} |
986 |
+ |
|
987 |
+ |
JunctionSet::const_iterator junc_iter = temp_junctions.begin(); |
988 |
+ |
for (; junc_iter != temp_junctions.end(); ++junc_iter) |
989 |
+ |
{ |
990 |
+ |
Junction junc = junc_iter->first; |
991 |
+ |
|
992 |
+ |
/* |
993 |
+ |
fprintf(stderr, "(type%d) %d %d-%d %s (AS:%d XM:%d) with junc %u-%u\n", |
994 |
+ |
!use_rev_junctions, |
995 |
+ |
bh.insert_id(), bh.left(), bh.right(), |
996 |
+ |
print_cigar(bh.cigar()).c_str(), |
997 |
+ |
bh.alignment_score(), bh.edit_dist(), |
998 |
+ |
junc.left, junc.right); |
999 |
+ |
//*/ |
1000 |
+ |
|
1001 |
+ |
int new_left = bh.left(); |
1002 |
+ |
int intron_length = junc.right - junc.left - 1; |
1003 |
+ |
vector<CigarOp> new_cigars; |
1004 |
+ |
bool anchored = false; |
1005 |
+ |
if (j == 0 && bh.left() > (int)junc.left) |
1006 |
+ |
{ |
1007 |
+ |
new_left -= intron_length; |
1008 |
+ |
int before_match_length = junc.left - new_left + 1;; |
1009 |
+ |
int after_match_length = op.length - before_match_length; |
1010 |
+ |
|
1011 |
+ |
if (before_match_length > 0 && after_match_length) |
1012 |
+ |
{ |
1013 |
+ |
anchored = true; |
1014 |
+ |
new_cigars.push_back(CigarOp(MATCH, before_match_length)); |
1015 |
+ |
new_cigars.push_back(CigarOp(REF_SKIP, intron_length)); |
1016 |
+ |
new_cigars.push_back(CigarOp(MATCH, after_match_length)); |
1017 |
+ |
|
1018 |
+ |
new_cigars.insert(new_cigars.end(), cigars.begin() + 1, cigars.end()); |
1019 |
+ |
} |
1020 |
+ |
} |
1021 |
+ |
else |
1022 |
+ |
{ |
1023 |
+ |
new_cigars.insert(new_cigars.end(), cigars.begin(), cigars.end() - 1); |
1024 |
+ |
|
1025 |
+ |
int before_match_length = junc.left - pos + 1; |
1026 |
+ |
int after_match_length = op.length - before_match_length; |
1027 |
+ |
|
1028 |
+ |
if (before_match_length > 0 && after_match_length > 0) |
1029 |
+ |
{ |
1030 |
+ |
anchored = true; |
1031 |
+ |
new_cigars.push_back(CigarOp(MATCH, before_match_length)); |
1032 |
+ |
new_cigars.push_back(CigarOp(REF_SKIP, intron_length)); |
1033 |
+ |
new_cigars.push_back(CigarOp(MATCH, after_match_length)); |
1034 |
+ |
} |
1035 |
+ |
} |
1036 |
+ |
|
1037 |
+ |
if (anchored) |
1038 |
+ |
{ |
1039 |
+ |
BowtieHit new_bh(bh.ref_id(), |
1040 |
+ |
bh.ref_id2(), |
1041 |
+ |
bh.insert_id(), |
1042 |
+ |
new_left, |
1043 |
+ |
new_cigars, |
1044 |
+ |
bh.antisense_align(), |
1045 |
+ |
junc.antisense, |
1046 |
+ |
0, /* edit_dist - needs to be recalculated */ |
1047 |
+ |
0, /* splice_mms - needs to be recalculated */ |
1048 |
+ |
false); |
1049 |
+ |
|
1050 |
+ |
new_bh.seq(bh.seq()); |
1051 |
+ |
new_bh.qual(bh.qual()); |
1052 |
+ |
|
1053 |
+ |
vector<string> aux_fields; |
1054 |
+ |
bowtie_sam_extra(new_bh, rt, aux_fields); |
1055 |
+ |
|
1056 |
+ |
vector<string>::const_iterator aux_iter = aux_fields.begin(); |
1057 |
+ |
for (; aux_iter != aux_fields.end(); ++aux_iter) |
1058 |
+ |
{ |
1059 |
+ |
const string& aux_field = *aux_iter; |
1060 |
+ |
if (strncmp(aux_field.c_str(), "AS", 2) == 0) |
1061 |
+ |
{ |
1062 |
+ |
int alignment_score = atoi(aux_field.c_str() + 5); |
1063 |
+ |
new_bh.alignment_score(alignment_score); |
1064 |
+ |
} |
1065 |
+ |
else if (strncmp(aux_field.c_str(), "XM", 2) == 0) |
1066 |
+ |
{ |
1067 |
+ |
int XM_value = atoi(aux_field.c_str() + 5); |
1068 |
+ |
new_bh.edit_dist(XM_value); |
1069 |
+ |
} |
1070 |
+ |
} |
1071 |
+ |
|
1072 |
+ |
vector<string> sam_toks; |
1073 |
+ |
tokenize(bh.hitfile_rec().c_str(), "\t", sam_toks); |
1074 |
+ |
|
1075 |
+ |
char coord[20] = {0,}; |
1076 |
+ |
sprintf(coord, "%d", new_bh.left() + 1); |
1077 |
+ |
sam_toks[3] = coord; |
1078 |
+ |
sam_toks[5] = print_cigar(new_bh.cigar()); |
1079 |
+ |
for (size_t a = 11; a < sam_toks.size(); ++a) |
1080 |
+ |
{ |
1081 |
+ |
string& sam_tok = sam_toks[a]; |
1082 |
+ |
for (size_t b = 0; b < aux_fields.size(); ++b) |
1083 |
+ |
{ |
1084 |
+ |
const string& aux_tok = aux_fields[b]; |
1085 |
+ |
if (strncmp(sam_tok.c_str(), aux_tok.c_str(), 5) == 0) |
1086 |
+ |
{ |
1087 |
+ |
sam_tok = aux_tok; |
1088 |
+ |
break; |
1089 |
+ |
} |
1090 |
+ |
} |
1091 |
+ |
} |
1092 |
+ |
|
1093 |
+ |
if (!bh.is_spliced()) |
1094 |
+ |
{ |
1095 |
+ |
if (junc.antisense) |
1096 |
+ |
sam_toks.push_back("XS:A:-"); |
1097 |
+ |
else |
1098 |
+ |
sam_toks.push_back("XS:A:+"); |
1099 |
+ |
} |
1100 |
+ |
|
1101 |
+ |
|
1102 |
+ |
string new_rec = ""; |
1103 |
+ |
for (size_t d = 0; d < sam_toks.size(); ++d) |
1104 |
+ |
{ |
1105 |
+ |
new_rec += sam_toks[d]; |
1106 |
+ |
if (d < sam_toks.size() - 1) |
1107 |
+ |
new_rec += "\t"; |
1108 |
+ |
} |
1109 |
+ |
|
1110 |
+ |
new_bh.hitfile_rec(new_rec); |
1111 |
+ |
|
1112 |
+ |
if (new_bh.edit_dist() <= bh.edit_dist()) |
1113 |
+ |
additional_hits.push_back(new_bh); |
1114 |
+ |
|
1115 |
+ |
/* |
1116 |
+ |
fprintf(stderr, "\t%d %d-%d %s (AS:%d XM:%d) with junc %u-%u\n", |
1117 |
+ |
new_bh.insert_id(), new_bh.left(), new_bh.right(), |
1118 |
+ |
print_cigar(new_bh.cigar()).c_str(), |
1119 |
+ |
new_bh.alignment_score(), new_bh.edit_dist(), |
1120 |
+ |
junc.left, junc.right); |
1121 |
+ |
//*/ |
1122 |
+ |
} |
1123 |
+ |
} |
1124 |
+ |
} |
1125 |
+ |
|
1126 |
+ |
#if 0 |
1127 |
+ |
{ |
1128 |
+ |
DeletionSet::const_iterator lb, ub; |
1129 |
+ |
bool use_rev_deletions = (j == 0); |
1130 |
+ |
const DeletionSet& curr_deletions = (use_rev_deletions ? rev_deletions : deletions); |
1131 |
+ |
if (use_rev_deletions) |
1132 |
+ |
{ |
1133 |
+ |
lb = curr_deletions.upper_bound(Deletion(refid, left1, 0, true)); |
1134 |
+ |
ub = curr_deletions.lower_bound(Deletion(refid, left2, left2, false)); |
1135 |
+ |
} |
1136 |
+ |
else |
1137 |
+ |
{ |
1138 |
+ |
int common_right = left2 + 100; |
1139 |
+ |
lb = curr_deletions.upper_bound(Deletion(refid, left1, common_right, true)); |
1140 |
+ |
ub = curr_deletions.lower_bound(Deletion(refid, left2, common_right, false)); |
1141 |
+ |
} |
1142 |
+ |
|
1143 |
+ |
while (lb != curr_deletions.end() && lb != ub) |
1144 |
+ |
{ |
1145 |
+ |
Deletion del = lb->first; |
1146 |
+ |
if (use_rev_deletions) |
1147 |
+ |
{ |
1148 |
+ |
int temp = del.left; |
1149 |
+ |
del.left = del.right; |
1150 |
+ |
del.right = temp; |
1151 |
+ |
} |
1152 |
+ |
|
1153 |
+ |
// daehwan - for debuggin purposes |
1154 |
+ |
/* |
1155 |
+ |
fprintf(stderr, "(type%d) %d %d-%d %s (AS:%d XM:%d) with junc %u-%u\n", |
1156 |
+ |
!use_rev_junctions, |
1157 |
+ |
bh.insert_id(), bh.left(), bh.right(), |
1158 |
+ |
print_cigar(bh.cigar()).c_str(), |
1159 |
+ |
bh.alignment_score(), bh.edit_dist(), |
1160 |
+ |
junc.left, junc.right); |
1161 |
+ |
*/ |
1162 |
+ |
|
1163 |
+ |
int del_length = del.right - del.left - 1; |
1164 |
+ |
int new_left = bh.left(); |
1165 |
+ |
if (j == 0) |
1166 |
+ |
new_left -= del_length; |
1167 |
+ |
|
1168 |
+ |
vector<CigarOp> new_cigars; |
1169 |
+ |
if (j == 0) |
1170 |
+ |
{ |
1171 |
+ |
int before_match_length = del.left - new_left + 1;; |
1172 |
+ |
int after_match_length = op.length - before_match_length; |
1173 |
+ |
|
1174 |
+ |
if (before_match_length > 0) |
1175 |
+ |
new_cigars.push_back(CigarOp(MATCH, before_match_length)); |
1176 |
+ |
new_cigars.push_back(CigarOp(DEL, del_length)); |
1177 |
+ |
if (after_match_length > 0) |
1178 |
+ |
new_cigars.push_back(CigarOp(MATCH, after_match_length)); |
1179 |
+ |
|
1180 |
+ |
new_cigars.insert(new_cigars.end(), cigars.begin() + 1, cigars.end()); |
1181 |
+ |
} |
1182 |
+ |
else |
1183 |
+ |
{ |
1184 |
+ |
new_cigars.insert(new_cigars.end(), cigars.begin(), cigars.end() - 1); |
1185 |
+ |
|
1186 |
+ |
int before_match_length = del.left - pos + 1; |
1187 |
+ |
int after_match_length = op.length - before_match_length; |
1188 |
+ |
|
1189 |
+ |
if (before_match_length > 0) |
1190 |
+ |
new_cigars.push_back(CigarOp(MATCH, before_match_length)); |
1191 |
+ |
new_cigars.push_back(CigarOp(DEL, del_length)); |
1192 |
+ |
if (after_match_length > 0) |
1193 |
+ |
new_cigars.push_back(CigarOp(MATCH, after_match_length)); |
1194 |
+ |
} |
1195 |
+ |
|
1196 |
+ |
BowtieHit new_bh(bh.ref_id(), |
1197 |
+ |
bh.ref_id2(), |
1198 |
+ |
bh.insert_id(), |
1199 |
+ |
new_left, |
1200 |
+ |
new_cigars, |
1201 |
+ |
bh.antisense_align(), |
1202 |
+ |
bh.antisense_splice(), |
1203 |
+ |
0, /* edit_dist - needs to be recalculated */ |
1204 |
+ |
0, /* splice_mms - needs to be recalculated */ |
1205 |
+ |
false); |
1206 |
+ |
|
1207 |
+ |
new_bh.seq(bh.seq()); |
1208 |
+ |
new_bh.qual(bh.qual()); |
1209 |
+ |
|
1210 |
+ |
vector<string> aux_fields; |
1211 |
+ |
bowtie_sam_extra(new_bh, rt, aux_fields); |
1212 |
+ |
|
1213 |
+ |
vector<string>::const_iterator aux_iter = aux_fields.begin(); |
1214 |
+ |
for (; aux_iter != aux_fields.end(); ++aux_iter) |
1215 |
+ |
{ |
1216 |
+ |
const string& aux_field = *aux_iter; |
1217 |
+ |
if (strncmp(aux_field.c_str(), "AS", 2) == 0) |
1218 |
+ |
{ |
1219 |
+ |
int alignment_score = atoi(aux_field.c_str() + 5); |
1220 |
+ |
new_bh.alignment_score(alignment_score); |
1221 |
+ |
} |
1222 |
+ |
else if (strncmp(aux_field.c_str(), "XM", 2) == 0) |
1223 |
+ |
{ |
1224 |
+ |
int XM_value = atoi(aux_field.c_str() + 5); |
1225 |
+ |
new_bh.edit_dist(XM_value); |
1226 |
+ |
} |
1227 |
+ |
} |
1228 |
+ |
|
1229 |
+ |
vector<string> sam_toks; |
1230 |
+ |
tokenize(bh.hitfile_rec().c_str(), "\t", sam_toks); |
1231 |
+ |
|
1232 |
+ |
char coord[20] = {0,}; |
1233 |
+ |
sprintf(coord, "%d", new_bh.left() + 1); |
1234 |
+ |
sam_toks[3] = coord; |
1235 |
+ |
sam_toks[5] = print_cigar(new_bh.cigar()); |
1236 |
+ |
for (size_t a = 11; a < sam_toks.size(); ++a) |
1237 |
+ |
{ |
1238 |
+ |
string& sam_tok = sam_toks[a]; |
1239 |
+ |
for (size_t b = 0; b < aux_fields.size(); ++b) |
1240 |
+ |
{ |
1241 |
+ |
const string& aux_tok = aux_fields[b]; |
1242 |
+ |
if (strncmp(sam_tok.c_str(), aux_tok.c_str(), 5) == 0) |
1243 |
+ |
{ |
1244 |
+ |
sam_tok = aux_tok; |
1245 |
+ |
break; |
1246 |
+ |
} |
1247 |
+ |
} |
1248 |
+ |
} |
1249 |
+ |
|
1250 |
+ |
string new_rec = ""; |
1251 |
+ |
for (size_t d = 0; d < sam_toks.size(); ++d) |
1252 |
+ |
{ |
1253 |
+ |
new_rec += sam_toks[d]; |
1254 |
+ |
if (d < sam_toks.size() - 1) |
1255 |
+ |
new_rec += "\t"; |
1256 |
+ |
} |
1257 |
+ |
|
1258 |
+ |
new_bh.hitfile_rec(new_rec); |
1259 |
+ |
|
1260 |
+ |
if (new_bh.edit_dist() <= bh.edit_dist()) |
1261 |
+ |
additional_hits.push_back(new_bh); |
1262 |
+ |
|
1263 |
+ |
/* |
1264 |
+ |
fprintf(stderr, "\t%d %d-%d %s (AS:%d XM:%d) with junc %u-%u\n", |
1265 |
+ |
new_bh.insert_id(), new_bh.left(), new_bh.right(), |
1266 |
+ |
print_cigar(new_bh.cigar()).c_str(), |
1267 |
+ |
new_bh.alignment_score(), new_bh.edit_dist(), |
1268 |
+ |
junc.left, junc.right); |
1269 |
+ |
*/ |
1270 |
+ |
|
1271 |
+ |
++lb; |
1272 |
+ |
} |
1273 |
+ |
} |
1274 |
+ |
|
1275 |
+ |
{ |
1276 |
+ |
InsertionSet::const_iterator lb, ub; |
1277 |
+ |
lb = insertions.upper_bound(Insertion(refid, left1, "")); |
1278 |
+ |
ub = insertions.lower_bound(Insertion(refid, left2, "")); |
1279 |
+ |
|
1280 |
+ |
while (lb != insertions.end() && lb != ub) |
1281 |
+ |
{ |
1282 |
+ |
// daehwan - for debugging purposse |
1283 |
+ |
break; |
1284 |
+ |
|
1285 |
+ |
Insertion ins = lb->first; |
1286 |
+ |
|
1287 |
+ |
// daehwan - for debugging purposes |
1288 |
+ |
/* |
1289 |
+ |
fprintf(stderr, "(type%d) %d %d-%d %s (AS:%d XM:%d) with junc %u-%u\n", |
1290 |
+ |
!use_rev_junctions, |
1291 |
+ |
bh.insert_id(), bh.left(), bh.right(), |
1292 |
+ |
print_cigar(bh.cigar()).c_str(), |
1293 |
+ |
bh.alignment_score(), bh.edit_dist(), |
1294 |
+ |
junc.left, junc.right); |
1295 |
+ |
*/ |
1296 |
+ |
|
1297 |
+ |
int ins_length = ins.sequence.length(); |
1298 |
+ |
int new_left = bh.left(); |
1299 |
+ |
if (j == 0) |
1300 |
+ |
new_left -= ins_length; |
1301 |
+ |
|
1302 |
+ |
vector<CigarOp> new_cigars; |
1303 |
+ |
if (j == 0) |
1304 |
+ |
{ |
1305 |
+ |
int before_match_length = ins.left - new_left + 1;; |
1306 |
+ |
int after_match_length = op.length - before_match_length - ins_length; |
1307 |
+ |
|
1308 |
+ |
if (before_match_length > 0) |
1309 |
+ |
new_cigars.push_back(CigarOp(MATCH, before_match_length)); |
1310 |
+ |
new_cigars.push_back(CigarOp(INS, ins_length)); |
1311 |
+ |
if (after_match_length > 0) |
1312 |
+ |
new_cigars.push_back(CigarOp(MATCH, after_match_length)); |
1313 |
+ |
|
1314 |
+ |
new_cigars.insert(new_cigars.end(), cigars.begin() + 1, cigars.end()); |
1315 |
+ |
} |
1316 |
+ |
else |
1317 |
+ |
{ |
1318 |
+ |
new_cigars.insert(new_cigars.end(), cigars.begin(), cigars.end() - 1); |
1319 |
+ |
|
1320 |
+ |
int before_match_length = ins.left - pos + 1; |
1321 |
+ |
int after_match_length = op.length - before_match_length - ins_length; |
1322 |
+ |
|
1323 |
+ |
if (before_match_length > 0) |
1324 |
+ |
new_cigars.push_back(CigarOp(MATCH, before_match_length)); |
1325 |
+ |
new_cigars.push_back(CigarOp(INS, ins_length)); |
1326 |
+ |
if (after_match_length > 0) |
1327 |
+ |
new_cigars.push_back(CigarOp(MATCH, after_match_length)); |
1328 |
+ |
} |
1329 |
+ |
|
1330 |
+ |
BowtieHit new_bh(bh.ref_id(), |
1331 |
+ |
bh.ref_id2(), |
1332 |
+ |
bh.insert_id(), |
1333 |
+ |
new_left, |
1334 |
+ |
new_cigars, |
1335 |
+ |
bh.antisense_align(), |
1336 |
+ |
bh.antisense_splice(), |
1337 |
+ |
0, /* edit_dist - needs to be recalculated */ |
1338 |
+ |
0, /* splice_mms - needs to be recalculated */ |
1339 |
+ |
false); |
1340 |
+ |
|
1341 |
+ |
new_bh.seq(bh.seq()); |
1342 |
+ |
new_bh.qual(bh.qual()); |
1343 |
+ |
|
1344 |
+ |
vector<string> aux_fields; |
1345 |
+ |
bowtie_sam_extra(new_bh, rt, aux_fields); |
1346 |
+ |
|
1347 |
+ |
vector<string>::const_iterator aux_iter = aux_fields.begin(); |
1348 |
+ |
for (; aux_iter != aux_fields.end(); ++aux_iter) |
1349 |
+ |
{ |
1350 |
+ |
const string& aux_field = *aux_iter; |
1351 |
+ |
if (strncmp(aux_field.c_str(), "AS", 2) == 0) |
1352 |
+ |
{ |
1353 |
+ |
int alignment_score = atoi(aux_field.c_str() + 5); |
1354 |
+ |
new_bh.alignment_score(alignment_score); |
1355 |
+ |
} |
1356 |
+ |
else if (strncmp(aux_field.c_str(), "XM", 2) == 0) |
1357 |
+ |
{ |
1358 |
+ |
int XM_value = atoi(aux_field.c_str() + 5); |
1359 |
+ |
new_bh.edit_dist(XM_value); |
1360 |
+ |
} |
1361 |
+ |
} |
1362 |
+ |
|
1363 |
+ |
/* |
1364 |
+ |
fprintf(stderr, "\t%d %d-%d %s (AS:%d XM:%d) with junc %u-%u\n", |
1365 |
+ |
new_bh.insert_id(), new_bh.left(), new_bh.right(), |
1366 |
+ |
print_cigar(new_bh.cigar()).c_str(), |
1367 |
+ |
new_bh.alignment_score(), new_bh.edit_dist(), |
1368 |
+ |
junc.left, junc.right); |
1369 |
+ |
*/ |
1370 |
+ |
|
1371 |
+ |
++lb; |
1372 |
+ |
} |
1373 |
+ |
} |
1374 |
+ |
#endif |
1375 |
+ |
} |
1376 |
+ |
} |
1377 |
+ |
|
1378 |
+ |
switch(op.opcode) |
1379 |
+ |
{ |
1380 |
+ |
case REF_SKIP: |
1381 |
+ |
pos += op.length; |
1382 |
+ |
break; |
1383 |
+ |
case rEF_SKIP: |
1384 |
+ |
pos -= op.length; |
1385 |
+ |
break; |
1386 |
+ |
case MATCH: |
1387 |
+ |
case DEL: |
1388 |
+ |
pos += op.length; |
1389 |
+ |
break; |
1390 |
+ |
case mATCH: |
1391 |
+ |
case dEL: |
1392 |
+ |
pos -= op.length; |
1393 |
+ |
break; |
1394 |
+ |
case FUSION_FF: |
1395 |
+ |
case FUSION_FR: |
1396 |
+ |
case FUSION_RF: |
1397 |
+ |
case FUSION_RR: |
1398 |
+ |
pos = op.length; |
1399 |
+ |
refid = bh.ref_id2(); |
1400 |
+ |
break; |
1401 |
+ |
default: |
1402 |
+ |
break; |
1403 |
+ |
} |
1404 |
+ |
} |
1405 |
+ |
} |
1406 |
+ |
|
1407 |
+ |
hits.hits.insert(hits.hits.end(), additional_hits.begin(), additional_hits.end()); |
1408 |
+ |
|
1409 |
+ |
std::sort(hits.hits.begin(), hits.hits.end()); |
1410 |
+ |
vector<BowtieHit>::iterator new_end = std::unique(hits.hits.begin(), hits.hits.end()); |
1411 |
+ |
hits.hits.erase(new_end, hits.hits.end()); |
1412 |
+ |
} |
1413 |
+ |
|
1414 |
+ |
|
1415 |
|
// events include splice junction, indels, and fusion points. |
1416 |
|
struct ConsensusEventsWorker |
1417 |
|
{ |
1556 |
|
{ |
1557 |
|
void operator()() |
1558 |
|
{ |
1559 |
+ |
rng.seed(1); |
1560 |
+ |
|
1561 |
|
ReadTable it; |
1562 |
|
GBamWriter bam_writer(bam_output_fname.c_str(), sam_header.c_str()); |
1563 |
|
|
1575 |
|
GBamWriter* left_um_out=new GBamWriter(left_um_fname.c_str(), sam_header.c_str()); |
1576 |
|
GBamWriter* right_um_out=NULL; |
1577 |
|
|
1062 |
– |
//if (left_um_out.openWrite(left_um_fname.c_str(), zpacker)==NULL) |
1063 |
– |
// err_die("Error: cannot open file %s for writing!\n",left_um_fname.c_str()); |
1578 |
|
ReadStream right_reads_file(right_reads_fname); |
1579 |
|
if (right_reads_offset > 0) |
1580 |
|
right_reads_file.seek(right_reads_offset); |
1627 |
|
{ |
1628 |
|
HitsForRead best_hits; |
1629 |
|
best_hits.insert_id = curr_left_obs_order; |
1630 |
< |
bool got_read=left_reads_file.getRead(curr_left_obs_order, l_read, |
1631 |
< |
reads_format, false, begin_id, end_id, |
1632 |
< |
left_um_out, false); |
1633 |
< |
//assert(got_read); |
1120 |
< |
|
1121 |
< |
if (right_reads_file.file()) { |
1122 |
< |
got_read=right_reads_file.getRead(curr_left_obs_order, r_read, |
1123 |
< |
reads_format, false, begin_id, end_id, |
1124 |
< |
right_um_out, true); |
1125 |
< |
//assert(got_read); |
1126 |
< |
} |
1127 |
< |
|
1630 |
> |
|
1631 |
> |
realign_reads(curr_left_hit_group, *rt, *junctions, *rev_junctions, |
1632 |
> |
*insertions, *deletions, *rev_deletions, *fusions); |
1633 |
> |
|
1634 |
|
exclude_hits_on_filtered_junctions(*junctions, curr_left_hit_group); |
1635 |
|
|
1636 |
|
// Process hits for left singleton, select best alignments |
1640 |
|
|
1641 |
|
if (best_hits.hits.size() > 0) |
1642 |
|
{ |
1643 |
< |
update_junctions(best_hits, *final_junctions); |
1644 |
< |
update_insertions_and_deletions(best_hits, *final_insertions, *final_deletions); |
1645 |
< |
update_fusions(best_hits, *rt, *final_fusions, *fusions); |
1646 |
< |
|
1647 |
< |
print_sam_for_single(*rt, |
1648 |
< |
best_hits, |
1649 |
< |
(right_map_fname.empty() ? FRAG_UNPAIRED : FRAG_LEFT), |
1650 |
< |
l_read, |
1651 |
< |
bam_writer, |
1652 |
< |
left_um_out); |
1643 |
> |
bool got_read = left_reads_file.getRead(curr_left_obs_order, l_read, |
1644 |
> |
reads_format, false, begin_id, end_id, |
1645 |
> |
left_um_out, false); |
1646 |
> |
assert (got_read); |
1647 |
> |
|
1648 |
> |
if (got_read) |
1649 |
> |
{ |
1650 |
> |
update_junctions(best_hits, *final_junctions); |
1651 |
> |
update_insertions_and_deletions(best_hits, *final_insertions, *final_deletions); |
1652 |
> |
update_fusions(best_hits, *rt, *final_fusions, *fusions); |
1653 |
> |
|
1654 |
> |
print_sam_for_single(*rt, |
1655 |
> |
best_hits, |
1656 |
> |
(right_map_fname.empty() ? FRAG_UNPAIRED : FRAG_LEFT), |
1657 |
> |
l_read.alt_name, |
1658 |
> |
bam_writer, |
1659 |
> |
rng); |
1660 |
> |
} |
1661 |
|
} |
1662 |
+ |
|
1663 |
|
// Get next hit group |
1664 |
|
left_hs.next_read_hits(curr_left_hit_group); |
1665 |
|
curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id); |
1675 |
|
|
1676 |
|
if (curr_right_obs_order >= begin_id) |
1677 |
|
{ |
1678 |
< |
bool got_read=right_reads_file.getRead(curr_right_obs_order, r_read, |
1679 |
< |
reads_format, false, begin_id, end_id, |
1165 |
< |
right_um_out, false); |
1166 |
< |
|
1167 |
< |
//assert(got_read); |
1168 |
< |
/* |
1169 |
< |
fprintf(right_um_out.file, "@%s #MAPPED#\n%s\n+\n%s\n", r_read.alt_name.c_str(), |
1170 |
< |
r_read.seq.c_str(), r_read.qual.c_str()); |
1171 |
< |
*/ |
1172 |
< |
got_read=left_reads_file.getRead(curr_right_obs_order, l_read, |
1173 |
< |
reads_format, false, begin_id, end_id, |
1174 |
< |
left_um_out, true); |
1175 |
< |
//assert(got_read); |
1176 |
< |
|
1678 |
> |
realign_reads(curr_right_hit_group, *rt, *junctions, *rev_junctions, |
1679 |
> |
*insertions, *deletions, *rev_deletions, *fusions); |
1680 |
|
exclude_hits_on_filtered_junctions(*junctions, curr_right_hit_group); |
1681 |
|
|
1682 |
|
// Process hit for right singleton, select best alignments |
1686 |
|
|
1687 |
|
if (best_hits.hits.size() > 0) |
1688 |
|
{ |
1689 |
< |
update_junctions(best_hits, *final_junctions); |
1690 |
< |
update_insertions_and_deletions(best_hits, *final_insertions, *final_deletions); |
1691 |
< |
update_fusions(best_hits, *rt, *final_fusions, *fusions); |
1692 |
< |
|
1693 |
< |
print_sam_for_single(*rt, |
1694 |
< |
best_hits, |
1695 |
< |
FRAG_RIGHT, |
1696 |
< |
r_read, |
1697 |
< |
bam_writer, |
1698 |
< |
right_um_out); |
1689 |
> |
bool got_read =right_reads_file.getRead(curr_right_obs_order, r_read, |
1690 |
> |
reads_format, false, begin_id, end_id, |
1691 |
> |
right_um_out, false); |
1692 |
> |
assert (got_read); |
1693 |
> |
|
1694 |
> |
if (got_read) |
1695 |
> |
{ |
1696 |
> |
update_junctions(best_hits, *final_junctions); |
1697 |
> |
update_insertions_and_deletions(best_hits, *final_insertions, *final_deletions); |
1698 |
> |
update_fusions(best_hits, *rt, *final_fusions, *fusions); |
1699 |
> |
|
1700 |
> |
print_sam_for_single(*rt, |
1701 |
> |
best_hits, |
1702 |
> |
FRAG_RIGHT, |
1703 |
> |
r_read.alt_name, |
1704 |
> |
bam_writer, |
1705 |
> |
rng); |
1706 |
> |
} |
1707 |
|
} |
1708 |
|
} |
1709 |
|
|
1718 |
|
curr_left_obs_order < end_id && |
1719 |
|
curr_left_obs_order != VMAXINT32) |
1720 |
|
{ |
1721 |
+ |
realign_reads(curr_left_hit_group, *rt, *junctions, *rev_junctions, |
1722 |
+ |
*insertions, *deletions, *rev_deletions, *fusions); |
1723 |
|
exclude_hits_on_filtered_junctions(*junctions, curr_left_hit_group); |
1724 |
+ |
|
1725 |
+ |
realign_reads(curr_right_hit_group, *rt, *junctions, *rev_junctions, |
1726 |
+ |
*insertions, *deletions, *rev_deletions, *fusions); |
1727 |
|
exclude_hits_on_filtered_junctions(*junctions, curr_right_hit_group); |
1728 |
+ |
|
1729 |
|
if (curr_left_hit_group.hits.empty()) |
1730 |
|
{ //only right read mapped |
1731 |
|
//write it in the mapped file with the #MAPPED# flag |
1215 |
– |
|
1216 |
– |
bool got_read=right_reads_file.getRead(curr_left_obs_order, r_read, |
1217 |
– |
reads_format, false, begin_id, end_id, |
1218 |
– |
right_um_out, false); |
1219 |
– |
//assert(got_read); |
1220 |
– |
/* |
1221 |
– |
fprintf(right_um_out.file, "@%s #MAPPED#\n%s\n+\n%s\n", r_read.alt_name.c_str(), |
1222 |
– |
r_read.seq.c_str(), r_read.qual.c_str()); |
1223 |
– |
*/ |
1732 |
|
HitsForRead right_best_hits; |
1733 |
|
right_best_hits.insert_id = curr_right_obs_order; |
1734 |
|
|
1738 |
|
|
1739 |
|
if (right_best_hits.hits.size() > 0) |
1740 |
|
{ |
1741 |
< |
update_junctions(right_best_hits, *final_junctions); |
1742 |
< |
update_insertions_and_deletions(right_best_hits, *final_insertions, *final_deletions); |
1743 |
< |
update_fusions(right_best_hits, *rt, *final_fusions, *fusions); |
1744 |
< |
|
1745 |
< |
print_sam_for_single(*rt, |
1746 |
< |
right_best_hits, |
1747 |
< |
FRAG_RIGHT, |
1748 |
< |
r_read, |
1749 |
< |
bam_writer, |
1750 |
< |
right_um_out); |
1741 |
> |
bool got_read = right_reads_file.getRead(curr_left_obs_order, r_read, |
1742 |
> |
reads_format, false, begin_id, end_id, |
1743 |
> |
right_um_out, false); |
1744 |
> |
assert (got_read); |
1745 |
> |
|
1746 |
> |
if (got_read) |
1747 |
> |
{ |
1748 |
> |
update_junctions(right_best_hits, *final_junctions); |
1749 |
> |
update_insertions_and_deletions(right_best_hits, *final_insertions, *final_deletions); |
1750 |
> |
update_fusions(right_best_hits, *rt, *final_fusions, *fusions); |
1751 |
> |
|
1752 |
> |
print_sam_for_single(*rt, |
1753 |
> |
right_best_hits, |
1754 |
> |
FRAG_RIGHT, |
1755 |
> |
r_read.alt_name, |
1756 |
> |
bam_writer, |
1757 |
> |
rng); |
1758 |
> |
} |
1759 |
|
} |
1760 |
|
} |
1761 |
|
else if (curr_right_hit_group.hits.empty()) |
1763 |
|
HitsForRead left_best_hits; |
1764 |
|
left_best_hits.insert_id = curr_left_obs_order; |
1765 |
|
//only left read mapped |
1250 |
– |
bool got_read=left_reads_file.getRead(curr_left_obs_order, l_read, |
1251 |
– |
reads_format, false, begin_id, end_id, |
1252 |
– |
left_um_out, false); |
1253 |
– |
|
1766 |
|
// Process hits for left singleton, select best alignments |
1767 |
|
read_best_alignments(curr_left_hit_group, left_best_hits, *gtf_junctions, |
1768 |
|
*junctions, *insertions, *deletions, *fusions, *coverage, |
1770 |
|
|
1771 |
|
if (left_best_hits.hits.size() > 0) |
1772 |
|
{ |
1773 |
< |
update_junctions(left_best_hits, *final_junctions); |
1774 |
< |
update_insertions_and_deletions(left_best_hits, *final_insertions, *final_deletions); |
1775 |
< |
update_fusions(left_best_hits, *rt, *final_fusions, *fusions); |
1776 |
< |
|
1777 |
< |
print_sam_for_single(*rt, |
1778 |
< |
left_best_hits, |
1779 |
< |
FRAG_LEFT, |
1780 |
< |
l_read, |
1781 |
< |
bam_writer, |
1782 |
< |
left_um_out); |
1773 |
> |
bool got_read =left_reads_file.getRead(curr_left_obs_order, l_read, |
1774 |
> |
reads_format, false, begin_id, end_id, |
1775 |
> |
left_um_out, false); |
1776 |
> |
assert (got_read); |
1777 |
> |
|
1778 |
> |
if (got_read) |
1779 |
> |
{ |
1780 |
> |
update_junctions(left_best_hits, *final_junctions); |
1781 |
> |
update_insertions_and_deletions(left_best_hits, *final_insertions, *final_deletions); |
1782 |
> |
update_fusions(left_best_hits, *rt, *final_fusions, *fusions); |
1783 |
> |
|
1784 |
> |
print_sam_for_single(*rt, |
1785 |
> |
left_best_hits, |
1786 |
> |
FRAG_LEFT, |
1787 |
> |
l_read.alt_name, |
1788 |
> |
bam_writer, |
1789 |
> |
rng); |
1790 |
> |
} |
1791 |
|
} |
1792 |
|
} |
1793 |
|
else |
1812 |
|
|
1813 |
|
if (best_hits.size() > 0) |
1814 |
|
{ |
1815 |
< |
update_junctions(best_left_hit_group, *final_junctions); |
1816 |
< |
update_insertions_and_deletions(best_left_hit_group, *final_insertions, *final_deletions); |
1817 |
< |
update_fusions(best_left_hit_group, *rt, *final_fusions, *fusions); |
1818 |
< |
|
1819 |
< |
update_junctions(best_right_hit_group, *final_junctions); |
1820 |
< |
update_insertions_and_deletions(best_right_hit_group, *final_insertions, *final_deletions); |
1821 |
< |
update_fusions(best_right_hit_group, *rt, *final_fusions, *fusions); |
1822 |
< |
|
1823 |
< |
pair_support(best_hits, *final_fusions, *fusions); |
1824 |
< |
|
1825 |
< |
print_sam_for_pair(*rt, |
1826 |
< |
best_hits, |
1827 |
< |
grade, |
1828 |
< |
left_reads_file, |
1829 |
< |
right_reads_file, |
1830 |
< |
bam_writer, |
1831 |
< |
left_um_out, |
1832 |
< |
right_um_out, |
1833 |
< |
begin_id, |
1834 |
< |
end_id); |
1815 |
> |
bool got_left_read = left_reads_file.getRead(best_hits[0].first.insert_id(), l_read, |
1816 |
> |
reads_format, false, begin_id, end_id, |
1817 |
> |
left_um_out, false); |
1818 |
> |
|
1819 |
> |
bool got_right_read = right_reads_file.getRead(best_hits[0].first.insert_id(), r_read, |
1820 |
> |
reads_format, false, begin_id, end_id, |
1821 |
> |
right_um_out, false); |
1822 |
> |
|
1823 |
> |
assert (got_left_read && got_right_read); |
1824 |
> |
|
1825 |
> |
if (got_left_read && got_right_read) |
1826 |
> |
{ |
1827 |
> |
update_junctions(best_left_hit_group, *final_junctions); |
1828 |
> |
update_insertions_and_deletions(best_left_hit_group, *final_insertions, *final_deletions); |
1829 |
> |
update_fusions(best_left_hit_group, *rt, *final_fusions, *fusions); |
1830 |
> |
|
1831 |
> |
update_junctions(best_right_hit_group, *final_junctions); |
1832 |
> |
update_insertions_and_deletions(best_right_hit_group, *final_insertions, *final_deletions); |
1833 |
> |
update_fusions(best_right_hit_group, *rt, *final_fusions, *fusions); |
1834 |
> |
|
1835 |
> |
pair_support(best_hits, *final_fusions, *fusions); |
1836 |
> |
|
1837 |
> |
print_sam_for_pair(*rt, |
1838 |
> |
best_hits, |
1839 |
> |
grade, |
1840 |
> |
bam_writer, |
1841 |
> |
l_read.alt_name, |
1842 |
> |
r_read.alt_name, |
1843 |
> |
rng, |
1844 |
> |
begin_id, |
1845 |
> |
end_id); |
1846 |
> |
} |
1847 |
|
} |
1848 |
|
} |
1849 |
|
|
1891 |
|
JunctionSet* gtf_junctions; |
1892 |
|
|
1893 |
|
JunctionSet* junctions; |
1894 |
+ |
JunctionSet* rev_junctions; |
1895 |
|
InsertionSet* insertions; |
1896 |
|
DeletionSet* deletions; |
1897 |
+ |
DeletionSet* rev_deletions; |
1898 |
|
FusionSet* fusions; |
1899 |
|
Coverage* coverage; |
1900 |
|
|
1911 |
|
int64_t left_map_offset; |
1912 |
|
int64_t right_reads_offset; |
1913 |
|
int64_t right_map_offset; |
1914 |
+ |
|
1915 |
+ |
boost::random::mt19937 rng; |
1916 |
|
}; |
1917 |
|
|
1918 |
|
void driver(const string& bam_output_fname, |
1930 |
|
num_threads = 1; |
1931 |
|
|
1932 |
|
RefSequenceTable rt(sam_header, true); |
1933 |
< |
|
1398 |
< |
if (fusion_search) |
1399 |
< |
get_seqs(ref_stream, rt, true); |
1933 |
> |
get_seqs(ref_stream, rt, true); |
1934 |
|
|
1935 |
|
srandom(1); |
1936 |
|
|
2047 |
|
DeletionSet& deletions = vdeletions[0]; |
2048 |
|
FusionSet& fusions = vfusions[0]; |
2049 |
|
Coverage& coverage = vcoverages[0]; |
2050 |
< |
for (size_t i = 1; i < num_threads; ++i) |
2050 |
> |
for (int i = 1; i < num_threads; ++i) |
2051 |
|
{ |
2052 |
|
merge_with(junctions, vjunctions[i]); |
2053 |
|
vjunctions[i].clear(); |
2066 |
|
} |
2067 |
|
|
2068 |
|
coverage.calculate_coverage(); |
2069 |
+ |
|
2070 |
+ |
JunctionSet rev_junctions; |
2071 |
+ |
JunctionSet::const_iterator junction_iter = junctions.begin(); |
2072 |
+ |
for (; junction_iter != junctions.end(); ++junction_iter) |
2073 |
+ |
{ |
2074 |
+ |
const Junction& junction = junction_iter->first; |
2075 |
+ |
Junction rev_junction = junction; |
2076 |
+ |
rev_junction.left = junction.right; |
2077 |
+ |
rev_junction.right = junction.left; |
2078 |
+ |
rev_junctions[rev_junction] = junction_iter->second; |
2079 |
+ |
} |
2080 |
+ |
|
2081 |
+ |
DeletionSet rev_deletions; |
2082 |
+ |
#if 0 |
2083 |
+ |
DeletionSet::const_iterator deletion_iter = deletions.begin(); |
2084 |
+ |
for (; deletion_iter != deletions.end(); ++deletion_iter) |
2085 |
+ |
{ |
2086 |
+ |
const Deletion& deletion = deletion_iter->first; |
2087 |
+ |
Deletion rev_deletion = deletion; |
2088 |
+ |
rev_deletion.left = deletion.right; |
2089 |
+ |
rev_deletion.right = deletion.left; |
2090 |
+ |
rev_deletions[rev_deletion] = deletion_iter->second; |
2091 |
+ |
} |
2092 |
+ |
#endif |
2093 |
|
|
2094 |
|
size_t num_unfiltered_juncs = junctions.size(); |
2095 |
|
fprintf(stderr, "Loaded %lu junctions\n", (long unsigned int) num_unfiltered_juncs); |
2148 |
|
|
2149 |
|
worker.gtf_junctions = >f_junctions; |
2150 |
|
worker.junctions = &junctions; |
2151 |
+ |
worker.rev_junctions = &rev_junctions; |
2152 |
|
worker.insertions = &insertions; |
2153 |
|
worker.deletions = &deletions; |
2154 |
+ |
worker.rev_deletions = &rev_deletions; |
2155 |
|
worker.fusions = &fusions; |
2156 |
|
worker.coverage = &coverage; |
2157 |
|
|
2204 |
|
InsertionSet final_insertions = vfinal_insertions[0]; |
2205 |
|
DeletionSet final_deletions = vfinal_deletions[0]; |
2206 |
|
FusionSet final_fusions = vfinal_fusions[0]; |
2207 |
< |
for (size_t i = 1; i < num_threads; ++i) |
2207 |
> |
for (int i = 1; i < num_threads; ++i) |
2208 |
|
{ |
2209 |
|
merge_with(final_junctions, vfinal_junctions[i]); |
2210 |
|
vfinal_junctions[i].clear(); |
2340 |
|
|
2341 |
|
string left_reads_filename = argv[optind++]; |
2342 |
|
string unzcmd=getUnpackCmd(left_reads_filename, false); |
1783 |
– |
|
2343 |
|
string right_map_filename; |
2344 |
|
string right_reads_filename; |
2345 |
|
|