| 1 |
/* |
| 2 |
* tophat_reports.cpp |
| 3 |
* TopHat |
| 4 |
* |
| 5 |
* Created by Cole Trapnell on 11/20/08. |
| 6 |
* Copyright 2008 Cole Trapnell. All rights reserved. |
| 7 |
* |
| 8 |
*/ |
| 9 |
|
| 10 |
#ifdef HAVE_CONFIG_H |
| 11 |
#include <config.h> |
| 12 |
#else |
| 13 |
#define PACKAGE_VERSION "INTERNAL" |
| 14 |
#define SVN_REVISION "XXX" |
| 15 |
#endif |
| 16 |
|
| 17 |
#include <cassert> |
| 18 |
#include <cstdio> |
| 19 |
#include <cstring> |
| 20 |
#include <vector> |
| 21 |
#include <string> |
| 22 |
#include <map> |
| 23 |
#include <algorithm> |
| 24 |
#include <set> |
| 25 |
#include <stdexcept> |
| 26 |
#include <iostream> |
| 27 |
#include <fstream> |
| 28 |
#include <seqan/sequence.h> |
| 29 |
#include <seqan/find.h> |
| 30 |
#include <seqan/file.h> |
| 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" |
| 38 |
#include "bwt_map.h" |
| 39 |
#include "junctions.h" |
| 40 |
#include "insertions.h" |
| 41 |
#include "deletions.h" |
| 42 |
#include "fusions.h" |
| 43 |
#include "align_status.h" |
| 44 |
#include "fragments.h" |
| 45 |
#include "wiggles.h" |
| 46 |
#include "tokenize.h" |
| 47 |
#include "reads.h" |
| 48 |
#include "coverage.h" |
| 49 |
|
| 50 |
|
| 51 |
#include "inserts.h" |
| 52 |
|
| 53 |
using namespace std; |
| 54 |
using namespace seqan; |
| 55 |
using std::set; |
| 56 |
|
| 57 |
static const JunctionSet empty_junctions; |
| 58 |
static const InsertionSet empty_insertions; |
| 59 |
static const DeletionSet empty_deletions; |
| 60 |
static const FusionSet empty_fusions; |
| 61 |
static const Coverage empty_coverage; |
| 62 |
|
| 63 |
// daehwan - this is redundancy, which should be removed. |
| 64 |
void get_seqs(istream& ref_stream, |
| 65 |
RefSequenceTable& rt, |
| 66 |
bool keep_seqs = true) |
| 67 |
{ |
| 68 |
while(ref_stream.good() && !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 |
seqan::read(ref_stream, *ref_str, Fasta()); |
| 79 |
|
| 80 |
rt.get_id(name, keep_seqs ? ref_str : NULL, 0); |
| 81 |
} |
| 82 |
} |
| 83 |
|
| 84 |
struct cmp_read_alignment |
| 85 |
{ |
| 86 |
bool operator() (const BowtieHit& l, const BowtieHit& r) const |
| 87 |
{ |
| 88 |
return l.alignment_score() > r.alignment_score(); |
| 89 |
} |
| 90 |
}; |
| 91 |
|
| 92 |
void read_best_alignments(const HitsForRead& hits_for_read, |
| 93 |
HitsForRead& best_hits, |
| 94 |
const JunctionSet& gtf_junctions, |
| 95 |
const JunctionSet& junctions = empty_junctions, |
| 96 |
const InsertionSet& insertions = empty_insertions, |
| 97 |
const DeletionSet& deletions = empty_deletions, |
| 98 |
const FusionSet& fusions = empty_fusions, |
| 99 |
const Coverage& coverage = empty_coverage, |
| 100 |
bool final_report = false) |
| 101 |
{ |
| 102 |
const vector<BowtieHit>& hits = hits_for_read.hits; |
| 103 |
|
| 104 |
if (hits.size() >= max_multihits * 5) |
| 105 |
return; |
| 106 |
|
| 107 |
for (size_t i = 0; i < hits.size(); ++i) |
| 108 |
{ |
| 109 |
if (hits[i].edit_dist() > max_read_mismatches) |
| 110 |
continue; |
| 111 |
|
| 112 |
BowtieHit hit = hits[i]; |
| 113 |
AlignStatus align_status(hit, gtf_junctions, |
| 114 |
junctions, insertions, deletions, fusions, coverage); |
| 115 |
hit.alignment_score(align_status._alignment_score); |
| 116 |
|
| 117 |
if (report_secondary_alignments || !final_report) |
| 118 |
{ |
| 119 |
best_hits.hits.push_back(hit); |
| 120 |
} |
| 121 |
else |
| 122 |
{ |
| 123 |
// Is the new status better than the current best one? |
| 124 |
if (best_hits.hits.size() == 0 || cmp_read_alignment()(hit, best_hits.hits[0])) |
| 125 |
{ |
| 126 |
best_hits.hits.clear(); |
| 127 |
best_hits.hits.push_back(hit); |
| 128 |
} |
| 129 |
else if (!cmp_read_alignment()(best_hits.hits[0], hit)) // is it just as good? |
| 130 |
{ |
| 131 |
best_hits.hits.push_back(hit); |
| 132 |
} |
| 133 |
} |
| 134 |
} |
| 135 |
|
| 136 |
if ((report_secondary_alignments || !final_report) && best_hits.hits.size() > 0) |
| 137 |
{ |
| 138 |
sort(best_hits.hits.begin(), best_hits.hits.end(), cmp_read_alignment()); |
| 139 |
} |
| 140 |
|
| 141 |
if (final_report) |
| 142 |
{ |
| 143 |
if (best_hits.hits.size() > max_multihits) |
| 144 |
best_hits.hits.erase(best_hits.hits.begin() + max_multihits, best_hits.hits.end()); |
| 145 |
} |
| 146 |
} |
| 147 |
|
| 148 |
bool set_insert_alignment_grade(const BowtieHit& lh, const BowtieHit& rh, InsertAlignmentGrade& grade) |
| 149 |
{ |
| 150 |
bool fusion = false; |
| 151 |
bool left_fusion = lh.fusion_opcode() != FUSION_NOTHING; |
| 152 |
bool right_fusion = rh.fusion_opcode() != FUSION_NOTHING; |
| 153 |
if (left_fusion && right_fusion) |
| 154 |
return false; |
| 155 |
|
| 156 |
fusion = left_fusion || right_fusion; |
| 157 |
if (!fusion && lh.ref_id() != rh.ref_id()) |
| 158 |
fusion = true; |
| 159 |
|
| 160 |
if (!fusion && lh.ref_id() == rh.ref_id()) |
| 161 |
{ |
| 162 |
if (lh.antisense_align() == rh.antisense_align()) |
| 163 |
fusion = true; |
| 164 |
else |
| 165 |
{ |
| 166 |
int inner_dist = 0; |
| 167 |
if (lh.antisense_align()) |
| 168 |
// used rh.left() instead of rh.right() for the cases, |
| 169 |
// where reads overlap with each other or reads span introns |
| 170 |
inner_dist = lh.left() - rh.left(); |
| 171 |
else |
| 172 |
inner_dist = rh.left() - lh.left(); |
| 173 |
|
| 174 |
if (inner_dist < 0 || inner_dist > (int)fusion_min_dist) |
| 175 |
fusion = true; |
| 176 |
} |
| 177 |
} |
| 178 |
|
| 179 |
grade = InsertAlignmentGrade(lh, rh, fusion); |
| 180 |
|
| 181 |
return true; |
| 182 |
} |
| 183 |
|
| 184 |
struct cmp_pair_alignment |
| 185 |
{ |
| 186 |
cmp_pair_alignment(const JunctionSet& gtf_juncs) : |
| 187 |
_gtf_juncs(>f_juncs) {} |
| 188 |
|
| 189 |
bool operator() (const pair<BowtieHit, BowtieHit>& l, const pair<BowtieHit, BowtieHit>& r) const |
| 190 |
{ |
| 191 |
InsertAlignmentGrade gl; set_insert_alignment_grade(l.first, l.second, gl); |
| 192 |
InsertAlignmentGrade gr; set_insert_alignment_grade(r.first, r.second, gr); |
| 193 |
|
| 194 |
bool better = gr < gl; |
| 195 |
bool worse = gl < gr; |
| 196 |
|
| 197 |
if (better && !worse) |
| 198 |
return true; |
| 199 |
else |
| 200 |
return false; |
| 201 |
} |
| 202 |
|
| 203 |
const JunctionSet* _gtf_juncs; |
| 204 |
}; |
| 205 |
|
| 206 |
void pair_best_alignments(const HitsForRead& left_hits, |
| 207 |
const HitsForRead& right_hits, |
| 208 |
InsertAlignmentGrade& best_grade, |
| 209 |
vector<pair<BowtieHit, BowtieHit> >& best_hits, |
| 210 |
const JunctionSet& gtf_junctions, |
| 211 |
const JunctionSet& junctions = empty_junctions, |
| 212 |
const InsertionSet& insertions = empty_insertions, |
| 213 |
const DeletionSet& deletions = empty_deletions, |
| 214 |
const FusionSet& fusions = empty_fusions, |
| 215 |
const Coverage& coverage = empty_coverage, |
| 216 |
bool final_report = false) |
| 217 |
{ |
| 218 |
const vector<BowtieHit>& left = left_hits.hits; |
| 219 |
const vector<BowtieHit>& right = right_hits.hits; |
| 220 |
|
| 221 |
if (left.size() >= max_multihits * 5 || right.size() >= max_multihits * 5) |
| 222 |
return; |
| 223 |
|
| 224 |
for (size_t i = 0; i < left.size(); ++i) |
| 225 |
{ |
| 226 |
if (left[i].edit_dist() > max_read_mismatches) continue; |
| 227 |
|
| 228 |
BowtieHit lh = left[i]; |
| 229 |
AlignStatus align_status(lh, gtf_junctions, |
| 230 |
junctions, insertions, deletions, fusions, coverage); |
| 231 |
lh.alignment_score(align_status._alignment_score); |
| 232 |
|
| 233 |
for (size_t j = 0; j < right.size(); ++j) |
| 234 |
{ |
| 235 |
if (right[j].edit_dist() > max_read_mismatches) continue; |
| 236 |
|
| 237 |
BowtieHit rh = right[j]; |
| 238 |
AlignStatus align_status(rh, gtf_junctions, |
| 239 |
junctions, insertions, deletions, fusions, coverage); |
| 240 |
rh.alignment_score(align_status._alignment_score); |
| 241 |
InsertAlignmentGrade g; |
| 242 |
bool allowed = set_insert_alignment_grade(lh, rh, g); |
| 243 |
|
| 244 |
// daehwan - for debugging purposes |
| 245 |
if (lh.insert_id() == 325708 && !g.is_fusion() && false) |
| 246 |
{ |
| 247 |
fprintf(stderr, "lh %d:%d %s score: %d (from %d) NM: %d\n", |
| 248 |
lh.ref_id(), lh.left(), print_cigar(lh.cigar()).c_str(), |
| 249 |
lh.alignment_score(), left[i].alignment_score(), lh.edit_dist()); |
| 250 |
fprintf(stderr, "rh %d:%d %s score: %d (from %d) NM: %d\n", |
| 251 |
rh.ref_id(), rh.left(), print_cigar(rh.cigar()).c_str(), |
| 252 |
rh.alignment_score(), right[j].alignment_score(), rh.edit_dist()); |
| 253 |
fprintf(stderr, "combined score: %d is_fusion(%d)\n", g.align_score(), g.is_fusion()); |
| 254 |
} |
| 255 |
|
| 256 |
if (!allowed) continue; |
| 257 |
if (!fusion_search && !report_discordant_pair_alignments && g.is_fusion()) continue; |
| 258 |
|
| 259 |
if (report_secondary_alignments || !final_report) |
| 260 |
{ |
| 261 |
best_hits.push_back(make_pair(lh, rh)); |
| 262 |
} |
| 263 |
else |
| 264 |
{ |
| 265 |
// Is the new status better than the current best one? |
| 266 |
if (best_grade < g) |
| 267 |
{ |
| 268 |
best_hits.clear(); |
| 269 |
best_grade = g; |
| 270 |
best_hits.push_back(make_pair(lh, rh)); |
| 271 |
} |
| 272 |
else if (!(g < best_grade)) |
| 273 |
{ |
| 274 |
best_hits.push_back(make_pair(lh, rh)); |
| 275 |
} |
| 276 |
} |
| 277 |
} |
| 278 |
} |
| 279 |
|
| 280 |
if ((report_secondary_alignments || !final_report) && best_hits.size() > 0) |
| 281 |
{ |
| 282 |
cmp_pair_alignment cmp(gtf_junctions); |
| 283 |
sort(best_hits.begin(), best_hits.end(), cmp); |
| 284 |
set_insert_alignment_grade(best_hits[0].first, best_hits[0].second, best_grade); |
| 285 |
} |
| 286 |
|
| 287 |
if (final_report) |
| 288 |
{ |
| 289 |
if (best_hits.size() > max_multihits) |
| 290 |
best_hits.erase(best_hits.begin() + max_multihits, best_hits.end()); |
| 291 |
} |
| 292 |
|
| 293 |
best_grade.num_alignments = best_hits.size(); |
| 294 |
} |
| 295 |
|
| 296 |
enum FragmentType {FRAG_UNPAIRED, FRAG_LEFT, FRAG_RIGHT}; |
| 297 |
|
| 298 |
void add_auxData(vector<string>& auxdata, vector<string>& sam_toks, |
| 299 |
const RefSequenceTable& rt,const BowtieHit& bh, FragmentType insert_side, |
| 300 |
int num_hits, const BowtieHit* next_hit, int hitIndex) { |
| 301 |
bool XS_found = false; |
| 302 |
if (sam_toks.size()>11) { |
| 303 |
|
| 304 |
for (size_t i=11;i<sam_toks.size();++i) { |
| 305 |
if (sam_toks[i].find("NH:i:")==string::npos && |
| 306 |
sam_toks[i].find("XS:i:")==string::npos) |
| 307 |
auxdata.push_back(sam_toks[i]); |
| 308 |
|
| 309 |
if (sam_toks[i].find("XS:A:")!=string::npos) |
| 310 |
XS_found = true; |
| 311 |
} |
| 312 |
} |
| 313 |
string aux("NH:i:"); |
| 314 |
str_appendInt(aux, num_hits); |
| 315 |
auxdata.push_back(aux); |
| 316 |
if (next_hit) { |
| 317 |
const char* nh_ref_name = "="; |
| 318 |
nh_ref_name = rt.get_name(next_hit->ref_id()); |
| 319 |
assert (nh_ref_name != NULL); |
| 320 |
bool same_contig=(next_hit->ref_id()==bh.ref_id()); |
| 321 |
aux="CC:Z:"; |
| 322 |
aux+= (same_contig ? "=" : nh_ref_name); |
| 323 |
auxdata.push_back(aux); |
| 324 |
aux="CP:i:"; |
| 325 |
int nh_gpos=next_hit->left() + 1; |
| 326 |
str_appendInt(aux, nh_gpos); |
| 327 |
auxdata.push_back(aux); |
| 328 |
} //has next_hit |
| 329 |
// FIXME: this code is still a bit brittle, because it contains no |
| 330 |
// consistency check that the mates are on opposite strands - a current protocol |
| 331 |
// requirement, and that the strand indicated by the alignment is consistent |
| 332 |
// with the orientation of the splices (though that should be handled upstream). |
| 333 |
if (!XS_found) { |
| 334 |
const string xs_f("XS:A:+"); |
| 335 |
const string xs_r("XS:A:-"); |
| 336 |
if (library_type == FR_FIRSTSTRAND) { |
| 337 |
if (insert_side == FRAG_LEFT || insert_side == FRAG_UNPAIRED) { |
| 338 |
if (bh.antisense_align()) |
| 339 |
auxdata.push_back(xs_f); |
| 340 |
else |
| 341 |
auxdata.push_back(xs_r); |
| 342 |
} |
| 343 |
else { |
| 344 |
if (bh.antisense_align()) |
| 345 |
auxdata.push_back(xs_r); |
| 346 |
else |
| 347 |
auxdata.push_back(xs_f); |
| 348 |
} |
| 349 |
} |
| 350 |
else if (library_type == FR_SECONDSTRAND) { |
| 351 |
if (insert_side == FRAG_LEFT || insert_side == FRAG_UNPAIRED){ |
| 352 |
if (bh.antisense_align()) |
| 353 |
auxdata.push_back(xs_r); |
| 354 |
else |
| 355 |
auxdata.push_back(xs_f); |
| 356 |
} |
| 357 |
else |
| 358 |
{ |
| 359 |
if (bh.antisense_align()) |
| 360 |
auxdata.push_back(xs_f); |
| 361 |
else |
| 362 |
auxdata.push_back(xs_r); |
| 363 |
} |
| 364 |
} |
| 365 |
} |
| 366 |
if (hitIndex >= 0) |
| 367 |
{ |
| 368 |
string aux("HI:i:"); |
| 369 |
str_appendInt(aux, hitIndex); |
| 370 |
auxdata.push_back(aux); |
| 371 |
} |
| 372 |
} |
| 373 |
|
| 374 |
bool rewrite_sam_record(GBamWriter& bam_writer, |
| 375 |
const RefSequenceTable& rt, |
| 376 |
const BowtieHit& bh, |
| 377 |
const char* bwt_buf, |
| 378 |
const char* read_alt_name, |
| 379 |
FragmentType insert_side, |
| 380 |
int num_hits, |
| 381 |
const BowtieHit* next_hit, |
| 382 |
bool primary, |
| 383 |
int hitIndex) |
| 384 |
{ |
| 385 |
// Rewrite this hit, filling in the alt name, mate mapping |
| 386 |
// and setting the pair flag |
| 387 |
vector<string> sam_toks; |
| 388 |
tokenize(bwt_buf, "\t", sam_toks); |
| 389 |
|
| 390 |
string ref_name = sam_toks[2], ref_name2 = ""; |
| 391 |
char cigar1[1024] = {0}, cigar2[1024] = {0}; |
| 392 |
string left_seq, right_seq, left_qual, right_qual; |
| 393 |
int left = -1, left1 = -1, left2 = -1; |
| 394 |
bool fusion_alignment = false; |
| 395 |
size_t XF_index = 0; |
| 396 |
for (size_t i = 11; i < sam_toks.size(); ++i) |
| 397 |
{ |
| 398 |
string& tok = sam_toks[i]; |
| 399 |
if (strncmp(tok.c_str(), "XF", 2) == 0) |
| 400 |
{ |
| 401 |
fusion_alignment = true; |
| 402 |
XF_index = i; |
| 403 |
|
| 404 |
vector<string> tuple_fields; |
| 405 |
tokenize(tok.c_str(), " ", tuple_fields); |
| 406 |
vector<string> contigs; |
| 407 |
tokenize(tuple_fields[1].c_str(), "-", contigs); |
| 408 |
if (contigs.size() >= 2) |
| 409 |
{ |
| 410 |
ref_name = contigs[0]; |
| 411 |
ref_name2 = contigs[1]; |
| 412 |
} |
| 413 |
|
| 414 |
extract_partial_hits(bh, tuple_fields[4].c_str(), tuple_fields[5].c_str(), |
| 415 |
cigar1, cigar2, left_seq, right_seq, |
| 416 |
left_qual, right_qual, left1, left2); |
| 417 |
|
| 418 |
break; |
| 419 |
} |
| 420 |
|
| 421 |
else if (strncmp(tok.c_str(), "AS", 2) == 0) |
| 422 |
{ |
| 423 |
char AS_score[128] = {0}; |
| 424 |
sprintf(AS_score, "AS:i:%d", bh.alignment_score()); |
| 425 |
tok = AS_score; |
| 426 |
} |
| 427 |
} |
| 428 |
|
| 429 |
string qname(read_alt_name); |
| 430 |
size_t slash_pos=qname.rfind('/'); |
| 431 |
if (slash_pos!=string::npos) |
| 432 |
qname.resize(slash_pos); |
| 433 |
//read_alt_name as QNAME |
| 434 |
int flag=atoi(sam_toks[1].c_str()); //FLAG |
| 435 |
if (insert_side != FRAG_UNPAIRED) { |
| 436 |
//flag = atoi(sam_toks[1].c_str()); |
| 437 |
// mark this as a singleton mate |
| 438 |
flag |= 0x0001; |
| 439 |
if (insert_side == FRAG_LEFT) |
| 440 |
flag |= 0x0040; |
| 441 |
else if (insert_side == FRAG_RIGHT) |
| 442 |
flag |= 0x0080; |
| 443 |
flag |= 0x0008; |
| 444 |
//char flag_buf[64]; |
| 445 |
//sprintf(flag_buf, "%d", flag); |
| 446 |
//sam_toks[t] = flag_buf; |
| 447 |
} |
| 448 |
if (!primary) |
| 449 |
flag |= 0x100; |
| 450 |
|
| 451 |
int gpos=isdigit(sam_toks[3][0]) ? atoi(sam_toks[3].c_str()) : 0; |
| 452 |
int mapQ = 50; |
| 453 |
if (num_hits > 1) { |
| 454 |
double err_prob = 1 - (1.0 / num_hits); |
| 455 |
mapQ = (int)(-10.0 * log(err_prob) / log(10.0)); |
| 456 |
} |
| 457 |
int tlen =atoi(sam_toks[8].c_str()); //TLEN |
| 458 |
int mate_pos=atoi(sam_toks[7].c_str()); |
| 459 |
|
| 460 |
GBamRecord* bamrec=NULL; |
| 461 |
if (fusion_alignment) { |
| 462 |
vector<string> auxdata; |
| 463 |
add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex); |
| 464 |
bamrec=bam_writer.new_record(qname.c_str(), flag, ref_name.c_str(), left1 + 1, mapQ, |
| 465 |
cigar1, sam_toks[6].c_str(), mate_pos, |
| 466 |
tlen, left_seq.c_str(), left_qual.c_str(), &auxdata); |
| 467 |
bam_writer.write(bamrec); |
| 468 |
delete bamrec; |
| 469 |
|
| 470 |
auxdata.clear(); |
| 471 |
sam_toks[XF_index][5] = '2'; |
| 472 |
add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex); |
| 473 |
bamrec=bam_writer.new_record(qname.c_str(), flag, ref_name2.c_str(), left2 + 1, mapQ, |
| 474 |
cigar2, sam_toks[6].c_str(), mate_pos, |
| 475 |
tlen, right_seq.c_str(), right_qual.c_str(), &auxdata); |
| 476 |
bam_writer.write(bamrec); |
| 477 |
delete bamrec; |
| 478 |
} else { |
| 479 |
vector<string> auxdata; |
| 480 |
add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex); |
| 481 |
bamrec=bam_writer.new_record(qname.c_str(), flag, sam_toks[2].c_str(), gpos, mapQ, |
| 482 |
sam_toks[5].c_str(), sam_toks[6].c_str(), mate_pos, |
| 483 |
tlen, sam_toks[9].c_str(), sam_toks[10].c_str(), &auxdata); |
| 484 |
bam_writer.write(bamrec); |
| 485 |
delete bamrec; |
| 486 |
} |
| 487 |
|
| 488 |
return true; |
| 489 |
} |
| 490 |
|
| 491 |
bool rewrite_sam_record(GBamWriter& bam_writer, |
| 492 |
const RefSequenceTable& rt, |
| 493 |
const BowtieHit& bh, |
| 494 |
const char* bwt_buf, |
| 495 |
const char* read_alt_name, |
| 496 |
const InsertAlignmentGrade& grade, |
| 497 |
FragmentType insert_side, |
| 498 |
const BowtieHit* partner, |
| 499 |
int num_hits, |
| 500 |
const BowtieHit* next_hit, |
| 501 |
bool primary, |
| 502 |
int hitIndex) |
| 503 |
{ |
| 504 |
// Rewrite this hit, filling in the alt name, mate mapping |
| 505 |
// and setting the pair flag |
| 506 |
vector<string> sam_toks; |
| 507 |
tokenize(bwt_buf, "\t", sam_toks); |
| 508 |
string qname(read_alt_name); |
| 509 |
size_t slash_pos=qname.rfind('/'); |
| 510 |
if (slash_pos!=string::npos) |
| 511 |
qname.resize(slash_pos); |
| 512 |
//read_alt_name as QNAME |
| 513 |
int flag = atoi(sam_toks[1].c_str()); |
| 514 |
// 0x0010 (strand of query) is assumed to be set correctly |
| 515 |
// to begin with |
| 516 |
flag |= 0x0001; //it must be paired |
| 517 |
if (insert_side == FRAG_LEFT) |
| 518 |
flag |= 0x0040; |
| 519 |
else if (insert_side == FRAG_RIGHT) |
| 520 |
flag |= 0x0080; |
| 521 |
if (!primary) |
| 522 |
flag |= 0x100; |
| 523 |
|
| 524 |
string ref_name = sam_toks[2], ref_name2 = ""; |
| 525 |
char cigar1[1024] = {0}, cigar2[1024] = {0}; |
| 526 |
string left_seq, right_seq, left_qual, right_qual; |
| 527 |
int left = -1, left1 = -1, left2 = -1; |
| 528 |
bool fusion_alignment = false; |
| 529 |
size_t XF_tok_idx = 11; |
| 530 |
for (; XF_tok_idx < sam_toks.size(); ++XF_tok_idx) |
| 531 |
{ |
| 532 |
string& tok = sam_toks[XF_tok_idx]; |
| 533 |
if (strncmp(tok.c_str(), "XF", 2) == 0) |
| 534 |
{ |
| 535 |
fusion_alignment = true; |
| 536 |
|
| 537 |
vector<string> tuple_fields; |
| 538 |
tokenize(tok.c_str(), " ", tuple_fields); |
| 539 |
vector<string> contigs; |
| 540 |
tokenize(tuple_fields[1].c_str(), "-", contigs); |
| 541 |
if (contigs.size() >= 2) |
| 542 |
{ |
| 543 |
ref_name = contigs[0]; |
| 544 |
ref_name2 = contigs[1]; |
| 545 |
} |
| 546 |
|
| 547 |
extract_partial_hits(bh, tuple_fields[4].c_str(), tuple_fields[5].c_str(), |
| 548 |
cigar1, cigar2, left_seq, right_seq, |
| 549 |
left_qual, right_qual, left1, left2); |
| 550 |
|
| 551 |
break; |
| 552 |
} |
| 553 |
|
| 554 |
else if (strncmp(tok.c_str(), "AS", 2) == 0) |
| 555 |
{ |
| 556 |
char AS_score[128] = {0}; |
| 557 |
sprintf(AS_score, "AS:i:%d", bh.alignment_score()); |
| 558 |
tok = AS_score; |
| 559 |
} |
| 560 |
} |
| 561 |
|
| 562 |
int gpos=isdigit(sam_toks[3][0]) ? atoi(sam_toks[3].c_str()) : 0; |
| 563 |
int mapQ = 50; |
| 564 |
if (grade.num_alignments > 1) { |
| 565 |
double err_prob = 1 - (1.0 / grade.num_alignments); |
| 566 |
mapQ = (int)(-10.0 * log(err_prob) / log(10.0)); |
| 567 |
} |
| 568 |
int tlen=0; //TLEN |
| 569 |
int mate_pos=atoi(sam_toks[7].c_str()); |
| 570 |
if (partner) { |
| 571 |
if (partner->ref_id()==bh.ref_id()) { |
| 572 |
sam_toks[6] = "="; //same chromosome |
| 573 |
//TLEN: |
| 574 |
tlen = bh.left() < partner->left() ? partner->right() - bh.left() : |
| 575 |
partner->left() - bh.right(); |
| 576 |
} |
| 577 |
|
| 578 |
else { //partner on different chromosome/contig |
| 579 |
sam_toks[6] = rt.get_name(partner->ref_id()); |
| 580 |
} |
| 581 |
mate_pos = partner->left() + 1; |
| 582 |
if (grade.happy()) |
| 583 |
flag |= 0x0002; |
| 584 |
if (partner->antisense_align()) |
| 585 |
flag |= 0x0020; |
| 586 |
|
| 587 |
if (partner->fusion_opcode() != FUSION_NOTHING) |
| 588 |
{ |
| 589 |
char partner_pos[1024]; |
| 590 |
sprintf(partner_pos, "XP:Z:%s-%s %d", rt.get_name(partner->ref_id()), rt.get_name(partner->ref_id2()), partner->left() + 1); |
| 591 |
sam_toks.push_back(partner_pos); |
| 592 |
} |
| 593 |
} |
| 594 |
else { |
| 595 |
sam_toks[6] = "*"; |
| 596 |
mate_pos = 0; |
| 597 |
flag |= 0x0008; |
| 598 |
} |
| 599 |
|
| 600 |
GBamRecord* bamrec=NULL; |
| 601 |
if (fusion_alignment) { |
| 602 |
vector<string> auxdata; |
| 603 |
add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex); |
| 604 |
bamrec=bam_writer.new_record(qname.c_str(), flag, ref_name.c_str(), left1 + 1, mapQ, |
| 605 |
cigar1, sam_toks[6].c_str(), mate_pos, |
| 606 |
tlen, left_seq.c_str(), left_qual.c_str(), &auxdata); |
| 607 |
bam_writer.write(bamrec); |
| 608 |
delete bamrec; |
| 609 |
|
| 610 |
auxdata.clear(); |
| 611 |
sam_toks[XF_tok_idx][5] = '2'; |
| 612 |
add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex); |
| 613 |
bamrec=bam_writer.new_record(qname.c_str(), flag, ref_name2.c_str(), left2 + 1, mapQ, |
| 614 |
cigar2, sam_toks[6].c_str(), mate_pos, |
| 615 |
tlen, right_seq.c_str(), right_qual.c_str(), &auxdata); |
| 616 |
bam_writer.write(bamrec); |
| 617 |
delete bamrec; |
| 618 |
} else { |
| 619 |
vector<string> auxdata; |
| 620 |
add_auxData(auxdata, sam_toks, rt, bh, insert_side, num_hits, next_hit, hitIndex); |
| 621 |
bamrec=bam_writer.new_record(qname.c_str(), flag, sam_toks[2].c_str(), gpos, mapQ, |
| 622 |
sam_toks[5].c_str(), sam_toks[6].c_str(), mate_pos, |
| 623 |
tlen, sam_toks[9].c_str(), sam_toks[10].c_str(), &auxdata); |
| 624 |
bam_writer.write(bamrec); |
| 625 |
delete bamrec; |
| 626 |
} |
| 627 |
|
| 628 |
return true; |
| 629 |
} |
| 630 |
|
| 631 |
struct lex_hit_sort |
| 632 |
{ |
| 633 |
lex_hit_sort(const RefSequenceTable& rt, const HitsForRead& hits) |
| 634 |
: _rt(rt), _hits(hits) |
| 635 |
{} |
| 636 |
|
| 637 |
bool operator()(const uint32_t& l, const uint32_t& r) const |
| 638 |
{ |
| 639 |
const BowtieHit& lhs = _hits.hits[l]; |
| 640 |
const BowtieHit& rhs = _hits.hits[r]; |
| 641 |
|
| 642 |
uint32_t l_id = lhs.ref_id(); |
| 643 |
uint32_t r_id = rhs.ref_id(); |
| 644 |
if (l_id != r_id) |
| 645 |
return l_id < r_id; |
| 646 |
|
| 647 |
return lhs.left() < rhs.left(); |
| 648 |
} |
| 649 |
|
| 650 |
const RefSequenceTable& _rt; |
| 651 |
const HitsForRead& _hits; |
| 652 |
}; |
| 653 |
|
| 654 |
void print_sam_for_single(const RefSequenceTable& rt, |
| 655 |
const HitsForRead& hits, |
| 656 |
FragmentType frag_type, |
| 657 |
Read& read, |
| 658 |
GBamWriter& bam_writer, |
| 659 |
GBamWriter* um_out, //write unmapped reads here, |
| 660 |
boost::random::mt19937& rng |
| 661 |
) |
| 662 |
{ |
| 663 |
assert(!read.alt_name.empty()); |
| 664 |
if (hits.hits.empty()) |
| 665 |
return; |
| 666 |
|
| 667 |
lex_hit_sort s(rt, hits); |
| 668 |
vector<uint32_t> index_vector; |
| 669 |
|
| 670 |
size_t i; |
| 671 |
for (i = 0; i < hits.hits.size(); ++i) |
| 672 |
index_vector.push_back(i); |
| 673 |
|
| 674 |
sort(index_vector.begin(), index_vector.end(), s); |
| 675 |
|
| 676 |
size_t primaryHit = 0; |
| 677 |
if (!report_secondary_alignments) |
| 678 |
primaryHit = rng() % hits.hits.size(); |
| 679 |
|
| 680 |
bool multipleHits = (hits.hits.size() > 1); |
| 681 |
for (i = 0; i < hits.hits.size(); ++i) |
| 682 |
{ |
| 683 |
size_t index = index_vector[i]; |
| 684 |
bool primary = (index == primaryHit); |
| 685 |
const BowtieHit& bh = hits.hits[index]; |
| 686 |
rewrite_sam_record(bam_writer, rt, |
| 687 |
bh, |
| 688 |
bh.hitfile_rec().c_str(), |
| 689 |
read.alt_name.c_str(), |
| 690 |
frag_type, |
| 691 |
hits.hits.size(), |
| 692 |
(i < hits.hits.size()-1) ? &(hits.hits[index_vector[i+1]]) : NULL, |
| 693 |
primary, |
| 694 |
(multipleHits? i: -1)); |
| 695 |
} |
| 696 |
} |
| 697 |
|
| 698 |
void print_sam_for_pair(const RefSequenceTable& rt, |
| 699 |
const vector<pair<BowtieHit, BowtieHit> >& best_hits, |
| 700 |
const InsertAlignmentGrade& grade, |
| 701 |
ReadStream& left_reads_file, |
| 702 |
ReadStream& right_reads_file, |
| 703 |
GBamWriter& bam_writer, |
| 704 |
GBamWriter* left_um_out, |
| 705 |
GBamWriter* right_um_out, |
| 706 |
boost::random::mt19937& rng, |
| 707 |
uint64_t begin_id = 0, |
| 708 |
uint64_t end_id = std::numeric_limits<uint64_t>::max()) |
| 709 |
{ |
| 710 |
Read left_read; |
| 711 |
Read right_read; |
| 712 |
if (best_hits.empty()) |
| 713 |
return; |
| 714 |
|
| 715 |
size_t i; |
| 716 |
HitsForRead right_hits; |
| 717 |
for (i = 0; i < best_hits.size(); ++i) |
| 718 |
right_hits.hits.push_back(best_hits[i].second); |
| 719 |
|
| 720 |
size_t primaryHit = 0; |
| 721 |
vector<uint32_t> index_vector; |
| 722 |
lex_hit_sort s(rt, right_hits); |
| 723 |
for (i = 0; i < right_hits.hits.size(); ++i) |
| 724 |
index_vector.push_back(i); |
| 725 |
|
| 726 |
sort(index_vector.begin(), index_vector.end(), s); |
| 727 |
|
| 728 |
if (!report_secondary_alignments) |
| 729 |
primaryHit = rng() % right_hits.hits.size(); |
| 730 |
|
| 731 |
bool got_left_read = left_reads_file.getRead(best_hits[0].first.insert_id(), left_read, |
| 732 |
reads_format, false, begin_id, end_id, |
| 733 |
left_um_out, false); |
| 734 |
|
| 735 |
bool got_right_read = right_reads_file.getRead(best_hits[0].first.insert_id(), right_read, |
| 736 |
reads_format, false, begin_id, end_id, |
| 737 |
right_um_out, false); |
| 738 |
|
| 739 |
assert (got_left_read && got_right_read); |
| 740 |
bool multipleHits = (best_hits.size() > 1); |
| 741 |
for (i = 0; i < best_hits.size(); ++i) |
| 742 |
{ |
| 743 |
size_t index = index_vector[i]; |
| 744 |
bool primary = (index == primaryHit); |
| 745 |
const BowtieHit& right_bh = best_hits[index].second; |
| 746 |
const BowtieHit& left_bh = best_hits[index].first; |
| 747 |
|
| 748 |
rewrite_sam_record(bam_writer, rt, |
| 749 |
right_bh, |
| 750 |
right_bh.hitfile_rec().c_str(), |
| 751 |
right_read.alt_name.c_str(), |
| 752 |
grade, |
| 753 |
FRAG_RIGHT, |
| 754 |
&left_bh, |
| 755 |
best_hits.size(), |
| 756 |
(i < best_hits.size() - 1) ? &(best_hits[index_vector[i+1]].second) : NULL, |
| 757 |
primary, |
| 758 |
(multipleHits? i: -1)); |
| 759 |
rewrite_sam_record(bam_writer, rt, |
| 760 |
left_bh, |
| 761 |
left_bh.hitfile_rec().c_str(), |
| 762 |
left_read.alt_name.c_str(), |
| 763 |
grade, |
| 764 |
FRAG_LEFT, |
| 765 |
&right_bh, |
| 766 |
best_hits.size(), |
| 767 |
(i < best_hits.size() - 1) ? &(best_hits[index_vector[i+1]].first) : NULL, |
| 768 |
primary, |
| 769 |
(multipleHits? i: -1)); |
| 770 |
} |
| 771 |
} |
| 772 |
|
| 773 |
/** |
| 774 |
* Given all of the hits for a particular read, update the read counts for insertions and deletions. |
| 775 |
* @param hits hits The alignments for a particular read |
| 776 |
* @param insertions Maps from an insertion to the number of supporting reads for that insertion |
| 777 |
* @param deletions Maps from a deletion to the number of supporting reads for that deletion |
| 778 |
*/ |
| 779 |
void update_insertions_and_deletions(const HitsForRead& hits, |
| 780 |
InsertionSet& insertions, |
| 781 |
DeletionSet& deletions) |
| 782 |
{ |
| 783 |
for (size_t i = 0; i < hits.hits.size(); ++i) |
| 784 |
{ |
| 785 |
const BowtieHit& bh = hits.hits[i]; |
| 786 |
insertions_from_alignment(bh, insertions); |
| 787 |
deletions_from_alignment(bh, deletions); |
| 788 |
} |
| 789 |
} |
| 790 |
|
| 791 |
void update_coverage(const HitsForRead& hits, |
| 792 |
Coverage& coverage) |
| 793 |
{ |
| 794 |
for (size_t i = 0; i < hits.hits.size(); ++i) |
| 795 |
{ |
| 796 |
const BowtieHit& hit = hits.hits[i]; |
| 797 |
const vector<CigarOp>& cigar = hit.cigar(); |
| 798 |
unsigned int positionInGenome = hit.left(); |
| 799 |
RefID ref_id = hit.ref_id(); |
| 800 |
|
| 801 |
for(size_t c = 0; c < cigar.size(); ++c) |
| 802 |
{ |
| 803 |
int opcode = cigar[c].opcode; |
| 804 |
int length = cigar[c].length; |
| 805 |
switch(opcode) |
| 806 |
{ |
| 807 |
case REF_SKIP: |
| 808 |
case MATCH: |
| 809 |
case DEL: |
| 810 |
if (opcode == MATCH) |
| 811 |
coverage.add_coverage(ref_id, positionInGenome, length); |
| 812 |
|
| 813 |
positionInGenome += length; |
| 814 |
break; |
| 815 |
case rEF_SKIP: |
| 816 |
case mATCH: |
| 817 |
case dEL: |
| 818 |
positionInGenome -= length; |
| 819 |
|
| 820 |
if (opcode == mATCH) |
| 821 |
coverage.add_coverage(ref_id, positionInGenome + 1, length); |
| 822 |
break; |
| 823 |
case FUSION_FF: |
| 824 |
case FUSION_FR: |
| 825 |
case FUSION_RF: |
| 826 |
case FUSION_RR: |
| 827 |
positionInGenome = length; |
| 828 |
ref_id = hit.ref_id2(); |
| 829 |
break; |
| 830 |
default: |
| 831 |
break; |
| 832 |
} |
| 833 |
} |
| 834 |
} |
| 835 |
} |
| 836 |
|
| 837 |
|
| 838 |
void update_fusions(const HitsForRead& hits, |
| 839 |
RefSequenceTable& rt, |
| 840 |
FusionSet& fusions, |
| 841 |
const FusionSet& fusions_ref = empty_fusions) |
| 842 |
{ |
| 843 |
if (!fusion_search) |
| 844 |
return; |
| 845 |
|
| 846 |
if (hits.hits.size() > fusion_multireads) |
| 847 |
return; |
| 848 |
|
| 849 |
bool update_stat = fusions_ref.size() > 0; |
| 850 |
for (size_t i = 0; i < hits.hits.size(); ++i) |
| 851 |
{ |
| 852 |
const BowtieHit& bh = hits.hits[i]; |
| 853 |
|
| 854 |
if (bh.edit_dist() > fusion_read_mismatches) |
| 855 |
continue; |
| 856 |
|
| 857 |
fusions_from_alignment(bh, fusions, rt, update_stat); |
| 858 |
|
| 859 |
if (update_stat) |
| 860 |
unsupport_fusions(bh, fusions, fusions_ref); |
| 861 |
} |
| 862 |
} |
| 863 |
|
| 864 |
void update_junctions(const HitsForRead& hits, |
| 865 |
JunctionSet& junctions) |
| 866 |
{ |
| 867 |
for (size_t i = 0; i < hits.hits.size(); ++i) |
| 868 |
{ |
| 869 |
const BowtieHit& bh = hits.hits[i]; |
| 870 |
junctions_from_alignment(bh, junctions); |
| 871 |
} |
| 872 |
} |
| 873 |
|
| 874 |
// Extracts junctions from all the SAM hits (based on REF_SKIPs) in the hit file |
| 875 |
// resets the stream when finished. |
| 876 |
void exclude_hits_on_filtered_junctions(const JunctionSet& junctions, HitsForRead& hits) |
| 877 |
{ |
| 878 |
HitsForRead remaining; |
| 879 |
remaining.insert_id = hits.insert_id; |
| 880 |
|
| 881 |
for (size_t i = 0; i < hits.hits.size(); ++i) |
| 882 |
{ |
| 883 |
BowtieHit& bh = hits.hits[i]; |
| 884 |
if (bh.edit_dist() > max_read_mismatches) |
| 885 |
continue; |
| 886 |
|
| 887 |
bool filter_hit = false; |
| 888 |
if (!bh.contiguous()) |
| 889 |
{ |
| 890 |
JunctionSet bh_juncs; |
| 891 |
junctions_from_alignment(bh, bh_juncs); |
| 892 |
for (JunctionSet::iterator itr = bh_juncs.begin(); |
| 893 |
itr != bh_juncs.end(); |
| 894 |
itr++) |
| 895 |
{ |
| 896 |
const Junction& j = itr->first; |
| 897 |
JunctionSet::const_iterator target = junctions.find(j); |
| 898 |
if (target == junctions.end() || !target->second.accepted) |
| 899 |
{ |
| 900 |
filter_hit = true; |
| 901 |
break; |
| 902 |
} |
| 903 |
} |
| 904 |
} |
| 905 |
if (!filter_hit) |
| 906 |
remaining.hits.push_back(bh); |
| 907 |
} |
| 908 |
hits = remaining; |
| 909 |
} |
| 910 |
|
| 911 |
// events include splice junction, indels, and fusion points. |
| 912 |
struct ConsensusEventsWorker |
| 913 |
{ |
| 914 |
void operator()() |
| 915 |
{ |
| 916 |
ReadTable it; |
| 917 |
vector<BAMHitFactory*> hit_factories; |
| 918 |
hit_factories.push_back(new BAMHitFactory(it, *rt)); |
| 919 |
HitStream l_hs(left_map_fname, hit_factories.back(), false, true, true, true); |
| 920 |
if (left_map_offset > 0) |
| 921 |
l_hs.seek(left_map_offset); |
| 922 |
|
| 923 |
hit_factories.push_back(new BAMHitFactory(it, *rt)); |
| 924 |
HitStream r_hs(right_map_fname, hit_factories.back(), false, true, true, true); |
| 925 |
if (right_map_offset > 0) |
| 926 |
r_hs.seek(right_map_offset); |
| 927 |
|
| 928 |
HitsForRead curr_left_hit_group; |
| 929 |
HitsForRead curr_right_hit_group; |
| 930 |
|
| 931 |
l_hs.next_read_hits(curr_left_hit_group); |
| 932 |
r_hs.next_read_hits(curr_right_hit_group); |
| 933 |
|
| 934 |
uint32_t curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id); |
| 935 |
uint32_t curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id); |
| 936 |
|
| 937 |
// While we still have unreported hits... |
| 938 |
while((curr_left_obs_order != VMAXINT32 || curr_right_obs_order != VMAXINT32) && |
| 939 |
(curr_left_obs_order < end_id || curr_right_obs_order < end_id)) |
| 940 |
{ |
| 941 |
// Chew up left singletons |
| 942 |
while (curr_left_obs_order < curr_right_obs_order && |
| 943 |
curr_left_obs_order < end_id && |
| 944 |
curr_left_obs_order != VMAXINT32) |
| 945 |
{ |
| 946 |
HitsForRead best_hits; |
| 947 |
best_hits.insert_id = curr_left_obs_order; |
| 948 |
|
| 949 |
// Process hits for left singleton, select best alignments |
| 950 |
read_best_alignments(curr_left_hit_group, best_hits, *gtf_junctions); |
| 951 |
update_coverage(best_hits, *coverage); |
| 952 |
update_junctions(best_hits, *junctions); |
| 953 |
update_insertions_and_deletions(best_hits, *insertions, *deletions); |
| 954 |
update_fusions(best_hits, *rt, *fusions); |
| 955 |
|
| 956 |
// Get next hit group |
| 957 |
l_hs.next_read_hits(curr_left_hit_group); |
| 958 |
curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id); |
| 959 |
} |
| 960 |
|
| 961 |
// Chew up right singletons |
| 962 |
while (curr_left_obs_order > curr_right_obs_order && |
| 963 |
curr_right_obs_order < end_id && |
| 964 |
curr_right_obs_order != VMAXINT32) |
| 965 |
{ |
| 966 |
HitsForRead best_hits; |
| 967 |
best_hits.insert_id = curr_right_obs_order; |
| 968 |
if (curr_right_obs_order >= begin_id) |
| 969 |
{ |
| 970 |
// Process hit for right singleton, select best alignments |
| 971 |
read_best_alignments(curr_right_hit_group, best_hits, *gtf_junctions); |
| 972 |
update_coverage(best_hits, *coverage); |
| 973 |
update_junctions(best_hits, *junctions); |
| 974 |
update_insertions_and_deletions(best_hits, *insertions, *deletions); |
| 975 |
update_fusions(best_hits, *rt, *fusions); |
| 976 |
} |
| 977 |
|
| 978 |
// Get next hit group |
| 979 |
r_hs.next_read_hits(curr_right_hit_group); |
| 980 |
curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id); |
| 981 |
} |
| 982 |
|
| 983 |
// Since we have both left hits and right hits for this insert, |
| 984 |
// Find the best pairing and print both |
| 985 |
while (curr_left_obs_order == curr_right_obs_order && |
| 986 |
curr_left_obs_order < end_id && |
| 987 |
curr_left_obs_order != VMAXINT32) |
| 988 |
{ |
| 989 |
vector<pair<BowtieHit, BowtieHit> > best_hits; |
| 990 |
|
| 991 |
InsertAlignmentGrade grade; |
| 992 |
pair_best_alignments(curr_left_hit_group, |
| 993 |
curr_right_hit_group, |
| 994 |
grade, |
| 995 |
best_hits, |
| 996 |
*gtf_junctions); |
| 997 |
|
| 998 |
HitsForRead best_left_hit_group; best_left_hit_group.insert_id = curr_left_obs_order; |
| 999 |
HitsForRead best_right_hit_group; best_right_hit_group.insert_id = curr_left_obs_order; |
| 1000 |
|
| 1001 |
for (size_t i = 0; i < best_hits.size(); ++i) |
| 1002 |
{ |
| 1003 |
best_left_hit_group.hits.push_back(best_hits[i].first); |
| 1004 |
best_right_hit_group.hits.push_back(best_hits[i].second); |
| 1005 |
} |
| 1006 |
|
| 1007 |
update_coverage(best_left_hit_group, *coverage); |
| 1008 |
update_junctions(best_left_hit_group, *junctions); |
| 1009 |
update_insertions_and_deletions(best_left_hit_group, *insertions, *deletions); |
| 1010 |
update_fusions(best_left_hit_group, *rt, *fusions); |
| 1011 |
|
| 1012 |
update_coverage(best_right_hit_group, *coverage); |
| 1013 |
update_junctions(best_right_hit_group, *junctions); |
| 1014 |
update_insertions_and_deletions(best_right_hit_group, *insertions, *deletions); |
| 1015 |
update_fusions(best_right_hit_group, *rt, *fusions); |
| 1016 |
|
| 1017 |
l_hs.next_read_hits(curr_left_hit_group); |
| 1018 |
curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id); |
| 1019 |
|
| 1020 |
r_hs.next_read_hits(curr_right_hit_group); |
| 1021 |
curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id); |
| 1022 |
} |
| 1023 |
} |
| 1024 |
|
| 1025 |
for (size_t i = 0; i < hit_factories.size(); ++i) |
| 1026 |
delete hit_factories[i]; |
| 1027 |
|
| 1028 |
hit_factories.clear(); |
| 1029 |
} |
| 1030 |
|
| 1031 |
string left_map_fname; |
| 1032 |
string right_map_fname; |
| 1033 |
|
| 1034 |
RefSequenceTable* rt; |
| 1035 |
|
| 1036 |
JunctionSet* gtf_junctions; |
| 1037 |
|
| 1038 |
uint64_t begin_id; |
| 1039 |
uint64_t end_id; |
| 1040 |
int64_t left_map_offset; |
| 1041 |
int64_t right_map_offset; |
| 1042 |
|
| 1043 |
JunctionSet* junctions; |
| 1044 |
InsertionSet* insertions; |
| 1045 |
DeletionSet* deletions; |
| 1046 |
FusionSet* fusions; |
| 1047 |
|
| 1048 |
Coverage* coverage; |
| 1049 |
}; |
| 1050 |
|
| 1051 |
struct ReportWorker |
| 1052 |
{ |
| 1053 |
void operator()() |
| 1054 |
{ |
| 1055 |
rng.seed(1); |
| 1056 |
|
| 1057 |
ReadTable it; |
| 1058 |
GBamWriter bam_writer(bam_output_fname.c_str(), sam_header.c_str()); |
| 1059 |
|
| 1060 |
ReadStream left_reads_file(left_reads_fname); |
| 1061 |
if (left_reads_file.file() == NULL) |
| 1062 |
err_die("Error: cannot open %s for reading\n", left_reads_fname.c_str()); |
| 1063 |
if (left_reads_file.isBam()) { |
| 1064 |
left_reads_file.use_alt_name(); |
| 1065 |
left_reads_file.ignoreQC(); |
| 1066 |
} |
| 1067 |
if (left_reads_offset > 0) |
| 1068 |
left_reads_file.seek(left_reads_offset); |
| 1069 |
|
| 1070 |
//if (!zpacker.empty()) left_um_fname+=".z"; |
| 1071 |
GBamWriter* left_um_out=new GBamWriter(left_um_fname.c_str(), sam_header.c_str()); |
| 1072 |
GBamWriter* right_um_out=NULL; |
| 1073 |
|
| 1074 |
ReadStream right_reads_file(right_reads_fname); |
| 1075 |
if (right_reads_offset > 0) |
| 1076 |
right_reads_file.seek(right_reads_offset); |
| 1077 |
|
| 1078 |
//FZPipe right_um_out; |
| 1079 |
if (!right_reads_fname.empty()) |
| 1080 |
{ |
| 1081 |
if (right_reads_file.isBam()) { |
| 1082 |
right_reads_file.use_alt_name(); |
| 1083 |
right_reads_file.ignoreQC(); |
| 1084 |
right_um_out=new GBamWriter(right_um_fname.c_str(), sam_header.c_str()); |
| 1085 |
} |
| 1086 |
//if (!zpacker.empty()) right_um_fname+=".z"; |
| 1087 |
//if (right_um_out.openWrite(right_um_fname.c_str(), zpacker)==NULL) |
| 1088 |
// err_die("Error: cannot open file %s for writing!\n",right_um_fname.c_str()); |
| 1089 |
} |
| 1090 |
|
| 1091 |
vector<BAMHitFactory*> hit_factories; |
| 1092 |
hit_factories.push_back(new BAMHitFactory(it, *rt)); |
| 1093 |
HitStream left_hs(left_map_fname, hit_factories.back(), false, true, true, true); |
| 1094 |
if (left_map_offset > 0) |
| 1095 |
left_hs.seek(left_map_offset); |
| 1096 |
|
| 1097 |
hit_factories.push_back(new BAMHitFactory(it, *rt)); |
| 1098 |
HitStream right_hs(right_map_fname, hit_factories.back(), false, true, true, true); |
| 1099 |
if (right_map_offset > 0) |
| 1100 |
right_hs.seek(right_map_offset); |
| 1101 |
|
| 1102 |
HitsForRead curr_left_hit_group; |
| 1103 |
HitsForRead curr_right_hit_group; |
| 1104 |
|
| 1105 |
left_hs.next_read_hits(curr_left_hit_group); |
| 1106 |
right_hs.next_read_hits(curr_right_hit_group); |
| 1107 |
|
| 1108 |
uint64_t curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id); |
| 1109 |
uint64_t curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id); |
| 1110 |
|
| 1111 |
const bool final_report = true; |
| 1112 |
|
| 1113 |
// While we still have unreported hits... |
| 1114 |
Read l_read; |
| 1115 |
Read r_read; |
| 1116 |
while((curr_left_obs_order != VMAXINT32 || curr_right_obs_order != VMAXINT32) && |
| 1117 |
(curr_left_obs_order < end_id || curr_right_obs_order < end_id)) |
| 1118 |
{ |
| 1119 |
// Chew up left singletons (pairs with right reads unmapped) |
| 1120 |
while (curr_left_obs_order < curr_right_obs_order && |
| 1121 |
curr_left_obs_order < end_id && |
| 1122 |
curr_left_obs_order != VMAXINT32) |
| 1123 |
{ |
| 1124 |
HitsForRead best_hits; |
| 1125 |
best_hits.insert_id = curr_left_obs_order; |
| 1126 |
bool got_read=left_reads_file.getRead(curr_left_obs_order, l_read, |
| 1127 |
reads_format, false, begin_id, end_id, |
| 1128 |
left_um_out, false); |
| 1129 |
//assert(got_read); |
| 1130 |
|
| 1131 |
if (right_reads_file.file()) { |
| 1132 |
got_read=right_reads_file.getRead(curr_left_obs_order, r_read, |
| 1133 |
reads_format, false, begin_id, end_id, |
| 1134 |
right_um_out, true); |
| 1135 |
//assert(got_read); |
| 1136 |
} |
| 1137 |
|
| 1138 |
exclude_hits_on_filtered_junctions(*junctions, curr_left_hit_group); |
| 1139 |
|
| 1140 |
// Process hits for left singleton, select best alignments |
| 1141 |
read_best_alignments(curr_left_hit_group, best_hits, *gtf_junctions, |
| 1142 |
*junctions, *insertions, *deletions, *fusions, *coverage, |
| 1143 |
final_report); |
| 1144 |
|
| 1145 |
if (best_hits.hits.size() > 0) |
| 1146 |
{ |
| 1147 |
update_junctions(best_hits, *final_junctions); |
| 1148 |
update_insertions_and_deletions(best_hits, *final_insertions, *final_deletions); |
| 1149 |
update_fusions(best_hits, *rt, *final_fusions, *fusions); |
| 1150 |
|
| 1151 |
print_sam_for_single(*rt, |
| 1152 |
best_hits, |
| 1153 |
(right_map_fname.empty() ? FRAG_UNPAIRED : FRAG_LEFT), |
| 1154 |
l_read, |
| 1155 |
bam_writer, |
| 1156 |
left_um_out, |
| 1157 |
rng); |
| 1158 |
} |
| 1159 |
else |
| 1160 |
{ |
| 1161 |
} |
| 1162 |
|
| 1163 |
// Get next hit group |
| 1164 |
left_hs.next_read_hits(curr_left_hit_group); |
| 1165 |
curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id); |
| 1166 |
} //left singletons |
| 1167 |
|
| 1168 |
// Chew up right singletons |
| 1169 |
while (curr_left_obs_order > curr_right_obs_order && |
| 1170 |
curr_right_obs_order < end_id && |
| 1171 |
curr_right_obs_order != VMAXINT32) |
| 1172 |
{ |
| 1173 |
HitsForRead best_hits; |
| 1174 |
best_hits.insert_id = curr_right_obs_order; |
| 1175 |
|
| 1176 |
if (curr_right_obs_order >= begin_id) |
| 1177 |
{ |
| 1178 |
bool got_read=right_reads_file.getRead(curr_right_obs_order, r_read, |
| 1179 |
reads_format, false, begin_id, end_id, |
| 1180 |
right_um_out, false); |
| 1181 |
|
| 1182 |
//assert(got_read); |
| 1183 |
/* |
| 1184 |
fprintf(right_um_out.file, "@%s #MAPPED#\n%s\n+\n%s\n", r_read.alt_name.c_str(), |
| 1185 |
r_read.seq.c_str(), r_read.qual.c_str()); |
| 1186 |
*/ |
| 1187 |
got_read=left_reads_file.getRead(curr_right_obs_order, l_read, |
| 1188 |
reads_format, false, begin_id, end_id, |
| 1189 |
left_um_out, true); |
| 1190 |
//assert(got_read); |
| 1191 |
|
| 1192 |
exclude_hits_on_filtered_junctions(*junctions, curr_right_hit_group); |
| 1193 |
|
| 1194 |
// Process hit for right singleton, select best alignments |
| 1195 |
read_best_alignments(curr_right_hit_group, best_hits, *gtf_junctions, |
| 1196 |
*junctions, *insertions, *deletions, *fusions, *coverage, |
| 1197 |
final_report); |
| 1198 |
|
| 1199 |
if (best_hits.hits.size() > 0) |
| 1200 |
{ |
| 1201 |
update_junctions(best_hits, *final_junctions); |
| 1202 |
update_insertions_and_deletions(best_hits, *final_insertions, *final_deletions); |
| 1203 |
update_fusions(best_hits, *rt, *final_fusions, *fusions); |
| 1204 |
|
| 1205 |
print_sam_for_single(*rt, |
| 1206 |
best_hits, |
| 1207 |
FRAG_RIGHT, |
| 1208 |
r_read, |
| 1209 |
bam_writer, |
| 1210 |
right_um_out, |
| 1211 |
rng); |
| 1212 |
} |
| 1213 |
else |
| 1214 |
{ |
| 1215 |
} |
| 1216 |
} |
| 1217 |
|
| 1218 |
// Get next hit group |
| 1219 |
right_hs.next_read_hits(curr_right_hit_group); |
| 1220 |
curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id); |
| 1221 |
} |
| 1222 |
|
| 1223 |
// Since we have both left hits and right hits for this insert, |
| 1224 |
// Find the best pairing and print both |
| 1225 |
while (curr_left_obs_order == curr_right_obs_order && |
| 1226 |
curr_left_obs_order < end_id && |
| 1227 |
curr_left_obs_order != VMAXINT32) |
| 1228 |
{ |
| 1229 |
exclude_hits_on_filtered_junctions(*junctions, curr_left_hit_group); |
| 1230 |
exclude_hits_on_filtered_junctions(*junctions, curr_right_hit_group); |
| 1231 |
if (curr_left_hit_group.hits.empty()) |
| 1232 |
{ //only right read mapped |
| 1233 |
//write it in the mapped file with the #MAPPED# flag |
| 1234 |
|
| 1235 |
bool got_read=right_reads_file.getRead(curr_left_obs_order, r_read, |
| 1236 |
reads_format, false, begin_id, end_id, |
| 1237 |
right_um_out, false); |
| 1238 |
//assert(got_read); |
| 1239 |
/* |
| 1240 |
fprintf(right_um_out.file, "@%s #MAPPED#\n%s\n+\n%s\n", r_read.alt_name.c_str(), |
| 1241 |
r_read.seq.c_str(), r_read.qual.c_str()); |
| 1242 |
*/ |
| 1243 |
HitsForRead right_best_hits; |
| 1244 |
right_best_hits.insert_id = curr_right_obs_order; |
| 1245 |
|
| 1246 |
read_best_alignments(curr_right_hit_group, right_best_hits, *gtf_junctions, |
| 1247 |
*junctions, *insertions, *deletions, *fusions, *coverage, |
| 1248 |
final_report); |
| 1249 |
|
| 1250 |
if (right_best_hits.hits.size() > 0) |
| 1251 |
{ |
| 1252 |
update_junctions(right_best_hits, *final_junctions); |
| 1253 |
update_insertions_and_deletions(right_best_hits, *final_insertions, *final_deletions); |
| 1254 |
update_fusions(right_best_hits, *rt, *final_fusions, *fusions); |
| 1255 |
|
| 1256 |
print_sam_for_single(*rt, |
| 1257 |
right_best_hits, |
| 1258 |
FRAG_RIGHT, |
| 1259 |
r_read, |
| 1260 |
bam_writer, |
| 1261 |
right_um_out, |
| 1262 |
rng); |
| 1263 |
} |
| 1264 |
} |
| 1265 |
else if (curr_right_hit_group.hits.empty()) |
| 1266 |
{ |
| 1267 |
HitsForRead left_best_hits; |
| 1268 |
left_best_hits.insert_id = curr_left_obs_order; |
| 1269 |
//only left read mapped |
| 1270 |
bool got_read=left_reads_file.getRead(curr_left_obs_order, l_read, |
| 1271 |
reads_format, false, begin_id, end_id, |
| 1272 |
left_um_out, false); |
| 1273 |
|
| 1274 |
// Process hits for left singleton, select best alignments |
| 1275 |
read_best_alignments(curr_left_hit_group, left_best_hits, *gtf_junctions, |
| 1276 |
*junctions, *insertions, *deletions, *fusions, *coverage, |
| 1277 |
final_report); |
| 1278 |
|
| 1279 |
if (left_best_hits.hits.size() > 0) |
| 1280 |
{ |
| 1281 |
update_junctions(left_best_hits, *final_junctions); |
| 1282 |
update_insertions_and_deletions(left_best_hits, *final_insertions, *final_deletions); |
| 1283 |
update_fusions(left_best_hits, *rt, *final_fusions, *fusions); |
| 1284 |
|
| 1285 |
print_sam_for_single(*rt, |
| 1286 |
left_best_hits, |
| 1287 |
FRAG_LEFT, |
| 1288 |
l_read, |
| 1289 |
bam_writer, |
| 1290 |
left_um_out, |
| 1291 |
rng); |
| 1292 |
} |
| 1293 |
} |
| 1294 |
else |
| 1295 |
{ |
| 1296 |
vector<pair<BowtieHit, BowtieHit> > best_hits; |
| 1297 |
|
| 1298 |
InsertAlignmentGrade grade; |
| 1299 |
pair_best_alignments(curr_left_hit_group, |
| 1300 |
curr_right_hit_group, |
| 1301 |
grade, best_hits, *gtf_junctions, |
| 1302 |
*junctions, *insertions, *deletions, *fusions, *coverage, |
| 1303 |
final_report); |
| 1304 |
|
| 1305 |
HitsForRead best_left_hit_group; best_left_hit_group.insert_id = curr_left_obs_order; |
| 1306 |
HitsForRead best_right_hit_group; best_right_hit_group.insert_id = curr_left_obs_order; |
| 1307 |
|
| 1308 |
for (size_t i = 0; i < best_hits.size(); ++i) |
| 1309 |
{ |
| 1310 |
best_left_hit_group.hits.push_back(best_hits[i].first); |
| 1311 |
best_right_hit_group.hits.push_back(best_hits[i].second); |
| 1312 |
} |
| 1313 |
|
| 1314 |
if (best_hits.size() > 0) |
| 1315 |
{ |
| 1316 |
update_junctions(best_left_hit_group, *final_junctions); |
| 1317 |
update_insertions_and_deletions(best_left_hit_group, *final_insertions, *final_deletions); |
| 1318 |
update_fusions(best_left_hit_group, *rt, *final_fusions, *fusions); |
| 1319 |
|
| 1320 |
update_junctions(best_right_hit_group, *final_junctions); |
| 1321 |
update_insertions_and_deletions(best_right_hit_group, *final_insertions, *final_deletions); |
| 1322 |
update_fusions(best_right_hit_group, *rt, *final_fusions, *fusions); |
| 1323 |
|
| 1324 |
pair_support(best_hits, *final_fusions, *fusions); |
| 1325 |
|
| 1326 |
print_sam_for_pair(*rt, |
| 1327 |
best_hits, |
| 1328 |
grade, |
| 1329 |
left_reads_file, |
| 1330 |
right_reads_file, |
| 1331 |
bam_writer, |
| 1332 |
left_um_out, |
| 1333 |
right_um_out, |
| 1334 |
rng, |
| 1335 |
begin_id, |
| 1336 |
end_id); |
| 1337 |
} |
| 1338 |
else |
| 1339 |
{ |
| 1340 |
} |
| 1341 |
} |
| 1342 |
|
| 1343 |
left_hs.next_read_hits(curr_left_hit_group); |
| 1344 |
curr_left_obs_order = it.observation_order(curr_left_hit_group.insert_id); |
| 1345 |
|
| 1346 |
right_hs.next_read_hits(curr_right_hit_group); |
| 1347 |
curr_right_obs_order = it.observation_order(curr_right_hit_group.insert_id); |
| 1348 |
} |
| 1349 |
} //while we still have unreported hits.. |
| 1350 |
//print the remaining unmapped reads at the end of each reads' stream |
| 1351 |
|
| 1352 |
left_reads_file.getRead(VMAXINT32, l_read, |
| 1353 |
reads_format, false, begin_id, end_id, |
| 1354 |
left_um_out, false); |
| 1355 |
if (right_reads_file.file()) |
| 1356 |
right_reads_file.getRead(VMAXINT32, r_read, |
| 1357 |
reads_format, false, begin_id, end_id, |
| 1358 |
right_um_out, false); |
| 1359 |
|
| 1360 |
|
| 1361 |
// pclose (pipe close), which waits for a process to end, seems to conflict with boost::thread::join somehow, |
| 1362 |
// resulting in deadlock like behavior. |
| 1363 |
delete left_um_out; |
| 1364 |
delete right_um_out; |
| 1365 |
|
| 1366 |
|
| 1367 |
for (size_t i = 0; i < hit_factories.size(); ++i) |
| 1368 |
delete hit_factories[i]; |
| 1369 |
|
| 1370 |
hit_factories.clear(); |
| 1371 |
} |
| 1372 |
|
| 1373 |
string bam_output_fname; |
| 1374 |
string sam_header_fname; |
| 1375 |
|
| 1376 |
string left_reads_fname; |
| 1377 |
string left_map_fname; |
| 1378 |
string right_reads_fname; |
| 1379 |
string right_map_fname; |
| 1380 |
|
| 1381 |
string left_um_fname; |
| 1382 |
string right_um_fname; |
| 1383 |
|
| 1384 |
JunctionSet* gtf_junctions; |
| 1385 |
|
| 1386 |
JunctionSet* junctions; |
| 1387 |
InsertionSet* insertions; |
| 1388 |
DeletionSet* deletions; |
| 1389 |
FusionSet* fusions; |
| 1390 |
Coverage* coverage; |
| 1391 |
|
| 1392 |
JunctionSet* final_junctions; |
| 1393 |
InsertionSet* final_insertions; |
| 1394 |
DeletionSet* final_deletions; |
| 1395 |
FusionSet* final_fusions; |
| 1396 |
|
| 1397 |
RefSequenceTable* rt; |
| 1398 |
|
| 1399 |
uint64_t begin_id; |
| 1400 |
uint64_t end_id; |
| 1401 |
int64_t left_reads_offset; |
| 1402 |
int64_t left_map_offset; |
| 1403 |
int64_t right_reads_offset; |
| 1404 |
int64_t right_map_offset; |
| 1405 |
|
| 1406 |
boost::random::mt19937 rng; |
| 1407 |
}; |
| 1408 |
|
| 1409 |
void driver(const string& bam_output_fname, |
| 1410 |
istream& ref_stream, |
| 1411 |
const string& left_map_fname, |
| 1412 |
const string& left_reads_fname, |
| 1413 |
const string& right_map_fname, |
| 1414 |
const string& right_reads_fname, |
| 1415 |
FILE* junctions_out, |
| 1416 |
FILE* insertions_out, |
| 1417 |
FILE* deletions_out, |
| 1418 |
FILE* fusions_out) |
| 1419 |
{ |
| 1420 |
if (!parallel) |
| 1421 |
num_threads = 1; |
| 1422 |
|
| 1423 |
RefSequenceTable rt(sam_header, true); |
| 1424 |
|
| 1425 |
if (fusion_search) |
| 1426 |
get_seqs(ref_stream, rt, true); |
| 1427 |
|
| 1428 |
srandom(1); |
| 1429 |
|
| 1430 |
JunctionSet gtf_junctions; |
| 1431 |
if (!gtf_juncs.empty()) |
| 1432 |
{ |
| 1433 |
char splice_buf[4096]; |
| 1434 |
FILE* splice_coords = fopen(gtf_juncs.c_str(), "r"); |
| 1435 |
if (splice_coords) |
| 1436 |
{ |
| 1437 |
while (fgets(splice_buf, sizeof(splice_buf), splice_coords)) |
| 1438 |
{ |
| 1439 |
char* nl = strrchr(splice_buf, '\n'); |
| 1440 |
char* buf = splice_buf; |
| 1441 |
if (nl) *nl = 0; |
| 1442 |
|
| 1443 |
char* ref_name = get_token((char**)&buf, "\t"); |
| 1444 |
char* scan_left_coord = get_token((char**)&buf, "\t"); |
| 1445 |
char* scan_right_coord = get_token((char**)&buf, "\t"); |
| 1446 |
char* orientation = get_token((char**)&buf, "\t"); |
| 1447 |
|
| 1448 |
if (!scan_left_coord || !scan_right_coord || !orientation) |
| 1449 |
{ |
| 1450 |
fprintf(stderr,"Error: malformed splice coordinate record in %s\n:%s\n", |
| 1451 |
gtf_juncs.c_str(), buf); |
| 1452 |
exit(1); |
| 1453 |
} |
| 1454 |
|
| 1455 |
uint32_t ref_id = rt.get_id(ref_name, NULL, 0); |
| 1456 |
uint32_t left_coord = atoi(scan_left_coord); |
| 1457 |
uint32_t right_coord = atoi(scan_right_coord); |
| 1458 |
bool antisense = *orientation == '-'; |
| 1459 |
|
| 1460 |
gtf_junctions.insert(make_pair<Junction, JunctionStats>(Junction(ref_id, left_coord, right_coord, antisense), JunctionStats())); |
| 1461 |
} |
| 1462 |
} |
| 1463 |
fprintf(stderr, "Loaded %d GFF junctions from %s.\n", (int)(gtf_junctions.size()), gtf_juncs.c_str()); |
| 1464 |
} |
| 1465 |
|
| 1466 |
vector<uint64_t> read_ids; |
| 1467 |
vector<vector<int64_t> > offsets; |
| 1468 |
if (num_threads > 1) |
| 1469 |
{ |
| 1470 |
vector<string> fnames; |
| 1471 |
if (right_map_fname != "") |
| 1472 |
{ |
| 1473 |
fnames.push_back(right_reads_fname); |
| 1474 |
fnames.push_back(right_map_fname); |
| 1475 |
} |
| 1476 |
fnames.push_back(left_reads_fname); |
| 1477 |
fnames.push_back(left_map_fname); |
| 1478 |
bool enough_data = calculate_offsets(fnames, read_ids, offsets); |
| 1479 |
if (!enough_data) |
| 1480 |
num_threads = 1; |
| 1481 |
} |
| 1482 |
|
| 1483 |
JunctionSet vjunctions[num_threads]; |
| 1484 |
InsertionSet vinsertions[num_threads]; |
| 1485 |
DeletionSet vdeletions[num_threads]; |
| 1486 |
FusionSet vfusions[num_threads]; |
| 1487 |
Coverage vcoverages[num_threads]; |
| 1488 |
|
| 1489 |
vector<boost::thread*> threads; |
| 1490 |
for (int i = 0; i < num_threads; ++i) |
| 1491 |
{ |
| 1492 |
ConsensusEventsWorker worker; |
| 1493 |
|
| 1494 |
worker.left_map_fname = left_map_fname; |
| 1495 |
worker.right_map_fname = right_map_fname; |
| 1496 |
worker.rt = &rt; |
| 1497 |
worker.gtf_junctions = >f_junctions; |
| 1498 |
|
| 1499 |
worker.junctions = &vjunctions[i]; |
| 1500 |
worker.insertions = &vinsertions[i]; |
| 1501 |
worker.deletions = &vdeletions[i]; |
| 1502 |
worker.fusions = &vfusions[i]; |
| 1503 |
worker.coverage = &vcoverages[i]; |
| 1504 |
|
| 1505 |
worker.right_map_offset = 0; |
| 1506 |
|
| 1507 |
if (i == 0) |
| 1508 |
{ |
| 1509 |
worker.begin_id = 0; |
| 1510 |
worker.left_map_offset = 0; |
| 1511 |
} |
| 1512 |
else |
| 1513 |
{ |
| 1514 |
size_t offsets_size = offsets[i-1].size(); |
| 1515 |
|
| 1516 |
worker.begin_id = read_ids[i-1]; |
| 1517 |
worker.left_map_offset = offsets[i-1].back(); |
| 1518 |
if (offsets_size == 4) |
| 1519 |
worker.right_map_offset = offsets[i-1][1]; |
| 1520 |
} |
| 1521 |
|
| 1522 |
worker.end_id = (i+1 < num_threads) ? read_ids[i] : std::numeric_limits<uint64_t>::max(); |
| 1523 |
|
| 1524 |
if (num_threads > 1) |
| 1525 |
threads.push_back(new boost::thread(worker)); |
| 1526 |
else |
| 1527 |
worker(); |
| 1528 |
} |
| 1529 |
|
| 1530 |
for (size_t i = 0; i < threads.size(); ++i) |
| 1531 |
{ |
| 1532 |
threads[i]->join(); |
| 1533 |
delete threads[i]; |
| 1534 |
threads[i] = NULL; |
| 1535 |
} |
| 1536 |
threads.clear(); |
| 1537 |
|
| 1538 |
JunctionSet& junctions = vjunctions[0]; |
| 1539 |
InsertionSet& insertions = vinsertions[0]; |
| 1540 |
DeletionSet& deletions = vdeletions[0]; |
| 1541 |
FusionSet& fusions = vfusions[0]; |
| 1542 |
Coverage& coverage = vcoverages[0]; |
| 1543 |
for (size_t i = 1; i < num_threads; ++i) |
| 1544 |
{ |
| 1545 |
merge_with(junctions, vjunctions[i]); |
| 1546 |
vjunctions[i].clear(); |
| 1547 |
|
| 1548 |
merge_with(insertions, vinsertions[i]); |
| 1549 |
vinsertions[i].clear(); |
| 1550 |
|
| 1551 |
merge_with(deletions, vdeletions[i]); |
| 1552 |
vdeletions[i].clear(); |
| 1553 |
|
| 1554 |
merge_with(fusions, vfusions[i]); |
| 1555 |
vfusions[i].clear(); |
| 1556 |
|
| 1557 |
coverage.merge_with(vcoverages[i]); |
| 1558 |
vcoverages[i].clear(); |
| 1559 |
} |
| 1560 |
|
| 1561 |
coverage.calculate_coverage(); |
| 1562 |
|
| 1563 |
size_t num_unfiltered_juncs = junctions.size(); |
| 1564 |
fprintf(stderr, "Loaded %lu junctions\n", (long unsigned int) num_unfiltered_juncs); |
| 1565 |
|
| 1566 |
// Read hits, extract junctions, and toss the ones that arent strongly enough supported. |
| 1567 |
filter_junctions(junctions, gtf_junctions); |
| 1568 |
|
| 1569 |
//size_t num_juncs_after_filter = junctions.size(); |
| 1570 |
//fprintf(stderr, "Filtered %lu junctions\n", |
| 1571 |
// num_unfiltered_juncs - num_juncs_after_filter); |
| 1572 |
|
| 1573 |
/* |
| 1574 |
size_t small_overhangs = 0; |
| 1575 |
for (JunctionSet::iterator i = junctions.begin(); i != junctions.end(); ++i) |
| 1576 |
{ |
| 1577 |
if (i->second.accepted && |
| 1578 |
(i->second.left_extent < min_anchor_len || i->second.right_extent < min_anchor_len)) |
| 1579 |
{ |
| 1580 |
small_overhangs++; |
| 1581 |
} |
| 1582 |
} |
| 1583 |
|
| 1584 |
if (small_overhangs >0) |
| 1585 |
fprintf(stderr, "Warning: %lu small overhang junctions!\n", (long unsigned int)small_overhangs); |
| 1586 |
*/ |
| 1587 |
|
| 1588 |
JunctionSet vfinal_junctions[num_threads]; |
| 1589 |
InsertionSet vfinal_insertions[num_threads]; |
| 1590 |
DeletionSet vfinal_deletions[num_threads]; |
| 1591 |
FusionSet vfinal_fusions[num_threads]; |
| 1592 |
|
| 1593 |
for (int i = 0; i < num_threads; ++i) |
| 1594 |
{ |
| 1595 |
ReportWorker worker; |
| 1596 |
|
| 1597 |
worker.sam_header_fname = sam_header; |
| 1598 |
char filename[1024] = {0}; |
| 1599 |
sprintf(filename, "%s%d.bam", bam_output_fname.c_str(), i); |
| 1600 |
worker.bam_output_fname = filename; |
| 1601 |
string tmpoutdir=getFdir(worker.bam_output_fname); |
| 1602 |
worker.left_um_fname = tmpoutdir; |
| 1603 |
sprintf(filename, "unmapped_left_%d.bam", i); |
| 1604 |
worker.left_um_fname+=filename; |
| 1605 |
|
| 1606 |
if (right_reads_fname != "") |
| 1607 |
{ |
| 1608 |
sprintf(filename, "unmapped_right_%d.bam", i); |
| 1609 |
worker.right_um_fname = tmpoutdir; |
| 1610 |
worker.right_um_fname += filename; |
| 1611 |
} |
| 1612 |
|
| 1613 |
worker.left_reads_fname = left_reads_fname; |
| 1614 |
worker.left_map_fname = left_map_fname; |
| 1615 |
worker.right_reads_fname = right_reads_fname; |
| 1616 |
worker.right_map_fname = right_map_fname; |
| 1617 |
|
| 1618 |
worker.gtf_junctions = >f_junctions; |
| 1619 |
worker.junctions = &junctions; |
| 1620 |
worker.insertions = &insertions; |
| 1621 |
worker.deletions = &deletions; |
| 1622 |
worker.fusions = &fusions; |
| 1623 |
worker.coverage = &coverage; |
| 1624 |
|
| 1625 |
worker.final_junctions = &vfinal_junctions[i]; |
| 1626 |
worker.final_insertions = &vfinal_insertions[i]; |
| 1627 |
worker.final_deletions = &vfinal_deletions[i]; |
| 1628 |
worker.final_fusions = &vfinal_fusions[i]; |
| 1629 |
worker.rt = &rt; |
| 1630 |
|
| 1631 |
worker.right_reads_offset = 0; |
| 1632 |
worker.right_map_offset = 0; |
| 1633 |
|
| 1634 |
if (i == 0) |
| 1635 |
{ |
| 1636 |
worker.begin_id = 0; |
| 1637 |
worker.left_reads_offset = 0; |
| 1638 |
worker.left_map_offset = 0; |
| 1639 |
} |
| 1640 |
else |
| 1641 |
{ |
| 1642 |
size_t offsets_size = offsets[i-1].size(); |
| 1643 |
|
| 1644 |
worker.begin_id = read_ids[i-1]; |
| 1645 |
worker.left_reads_offset = offsets[i-1][offsets_size - 2]; |
| 1646 |
worker.left_map_offset = offsets[i-1].back(); |
| 1647 |
if (offsets_size == 4) |
| 1648 |
{ |
| 1649 |
worker.right_reads_offset = offsets[i-1][0]; |
| 1650 |
worker.right_map_offset = offsets[i-1][1]; |
| 1651 |
} |
| 1652 |
} |
| 1653 |
|
| 1654 |
worker.end_id = (i+1 < num_threads) ? read_ids[i] : std::numeric_limits<uint64_t>::max(); |
| 1655 |
|
| 1656 |
if (num_threads > 1) |
| 1657 |
threads.push_back(new boost::thread(worker)); |
| 1658 |
else |
| 1659 |
worker(); |
| 1660 |
} |
| 1661 |
|
| 1662 |
for (size_t i = 0; i < threads.size(); ++i) |
| 1663 |
{ |
| 1664 |
threads[i]->join(); |
| 1665 |
delete threads[i]; |
| 1666 |
threads[i] = NULL; |
| 1667 |
} |
| 1668 |
threads.clear(); |
| 1669 |
|
| 1670 |
JunctionSet final_junctions = vfinal_junctions[0]; |
| 1671 |
InsertionSet final_insertions = vfinal_insertions[0]; |
| 1672 |
DeletionSet final_deletions = vfinal_deletions[0]; |
| 1673 |
FusionSet final_fusions = vfinal_fusions[0]; |
| 1674 |
for (size_t i = 1; i < num_threads; ++i) |
| 1675 |
{ |
| 1676 |
merge_with(final_junctions, vfinal_junctions[i]); |
| 1677 |
vfinal_junctions[i].clear(); |
| 1678 |
|
| 1679 |
merge_with(final_insertions, vfinal_insertions[i]); |
| 1680 |
vfinal_insertions[i].clear(); |
| 1681 |
|
| 1682 |
merge_with(final_deletions, vfinal_deletions[i]); |
| 1683 |
vfinal_deletions[i].clear(); |
| 1684 |
|
| 1685 |
merge_with(final_fusions, vfinal_fusions[i]); |
| 1686 |
vfinal_fusions[i].clear(); |
| 1687 |
} |
| 1688 |
|
| 1689 |
fprintf (stderr, "Reporting final accepted alignments..."); |
| 1690 |
fprintf (stderr, "done.\n"); |
| 1691 |
|
| 1692 |
//small_overhangs = 0; |
| 1693 |
for (JunctionSet::iterator i = final_junctions.begin(); i != final_junctions.end();) |
| 1694 |
{ |
| 1695 |
if (i->second.supporting_hits == 0 || i->second.left_extent < 8 || i->second.right_extent < 8) |
| 1696 |
{ |
| 1697 |
final_junctions.erase(i++); |
| 1698 |
} |
| 1699 |
else |
| 1700 |
{ |
| 1701 |
++i; |
| 1702 |
} |
| 1703 |
} |
| 1704 |
|
| 1705 |
// if (small_overhangs > 0) |
| 1706 |
// fprintf(stderr, "Warning: %lu small overhang junctions!\n", small_overhangs); |
| 1707 |
|
| 1708 |
fprintf (stderr, "Printing junction BED track..."); |
| 1709 |
print_junctions(junctions_out, final_junctions, rt); |
| 1710 |
fprintf (stderr, "done\n"); |
| 1711 |
|
| 1712 |
fprintf (stderr, "Printing insertions..."); |
| 1713 |
print_insertions(insertions_out, final_insertions,rt); |
| 1714 |
fclose(insertions_out); |
| 1715 |
fprintf (stderr, "done\n"); |
| 1716 |
|
| 1717 |
fprintf (stderr, "Printing deletions..."); |
| 1718 |
print_deletions(deletions_out, final_deletions, rt); |
| 1719 |
fclose(deletions_out); |
| 1720 |
fprintf (stderr, "done\n"); |
| 1721 |
|
| 1722 |
if (fusion_search) |
| 1723 |
{ |
| 1724 |
fprintf (stderr, "Printing fusions..."); |
| 1725 |
print_fusions(fusions_out, final_fusions, rt); |
| 1726 |
fclose(fusions_out); |
| 1727 |
fprintf (stderr, "done\n"); |
| 1728 |
} |
| 1729 |
|
| 1730 |
fprintf(stderr, "Found %lu junctions from happy spliced reads\n", (long unsigned int)final_junctions.size()); |
| 1731 |
} |
| 1732 |
|
| 1733 |
void print_usage() |
| 1734 |
{ |
| 1735 |
fprintf(stderr, "Usage: tophat_reports <junctions.bed> <insertions.vcf> <deletions.vcf> <accepted_hits.sam> <left_map1,...,left_mapN> <left_reads.fq> [right_map1,...,right_mapN] [right_reads.fq]\n"); |
| 1736 |
|
| 1737 |
// fprintf(stderr, "Usage: tophat_reports <coverage.wig> <junctions.bed> <accepted_hits.sam> <map1.bwtout> [splice_map1.sbwtout]\n"); |
| 1738 |
} |
| 1739 |
|
| 1740 |
int main(int argc, char** argv) |
| 1741 |
{ |
| 1742 |
fprintf(stderr, "tophat_reports v%s (%s)\n", PACKAGE_VERSION, SVN_REVISION); |
| 1743 |
fprintf(stderr, "---------------------------------------\n"); |
| 1744 |
|
| 1745 |
reads_format = FASTQ; |
| 1746 |
|
| 1747 |
int parse_ret = parse_options(argc, argv, print_usage); |
| 1748 |
if (parse_ret) |
| 1749 |
return parse_ret; |
| 1750 |
|
| 1751 |
if(optind >= argc) |
| 1752 |
{ |
| 1753 |
print_usage(); |
| 1754 |
return 1; |
| 1755 |
} |
| 1756 |
|
| 1757 |
string ref_file_name = argv[optind++]; |
| 1758 |
|
| 1759 |
if(optind >= argc) |
| 1760 |
{ |
| 1761 |
print_usage(); |
| 1762 |
return 1; |
| 1763 |
} |
| 1764 |
|
| 1765 |
string junctions_file_name = argv[optind++]; |
| 1766 |
|
| 1767 |
if(optind >= argc) |
| 1768 |
{ |
| 1769 |
print_usage(); |
| 1770 |
return 1; |
| 1771 |
} |
| 1772 |
|
| 1773 |
string insertions_file_name = argv[optind++]; |
| 1774 |
if(optind >= argc) |
| 1775 |
{ |
| 1776 |
print_usage(); |
| 1777 |
return 1; |
| 1778 |
} |
| 1779 |
|
| 1780 |
string deletions_file_name = argv[optind++]; |
| 1781 |
if(optind >= argc) |
| 1782 |
{ |
| 1783 |
print_usage(); |
| 1784 |
return 1; |
| 1785 |
} |
| 1786 |
|
| 1787 |
string fusions_file_name = argv[optind++]; |
| 1788 |
if(optind >= argc) |
| 1789 |
{ |
| 1790 |
print_usage(); |
| 1791 |
return 1; |
| 1792 |
} |
| 1793 |
|
| 1794 |
string accepted_hits_file_name = argv[optind++]; |
| 1795 |
if(optind >= argc) |
| 1796 |
{ |
| 1797 |
print_usage(); |
| 1798 |
return 1; |
| 1799 |
} |
| 1800 |
|
| 1801 |
string left_map_filename = argv[optind++]; |
| 1802 |
if(optind >= argc) |
| 1803 |
{ |
| 1804 |
print_usage(); |
| 1805 |
return 1; |
| 1806 |
} |
| 1807 |
|
| 1808 |
string left_reads_filename = argv[optind++]; |
| 1809 |
string unzcmd=getUnpackCmd(left_reads_filename, false); |
| 1810 |
string right_map_filename; |
| 1811 |
string right_reads_filename; |
| 1812 |
|
| 1813 |
if (optind < argc) |
| 1814 |
{ |
| 1815 |
right_map_filename = argv[optind++]; |
| 1816 |
if(optind >= argc) { |
| 1817 |
print_usage(); |
| 1818 |
return 1; |
| 1819 |
} |
| 1820 |
right_reads_filename=argv[optind++]; |
| 1821 |
} |
| 1822 |
|
| 1823 |
ifstream ref_stream(ref_file_name.c_str(), ifstream::in); |
| 1824 |
if (!ref_stream.good()) |
| 1825 |
{ |
| 1826 |
fprintf(stderr, "Error: cannot open %s for reading\n", |
| 1827 |
ref_file_name.c_str()); |
| 1828 |
exit(1); |
| 1829 |
} |
| 1830 |
|
| 1831 |
FILE* junctions_file = fopen(junctions_file_name.c_str(), "w"); |
| 1832 |
if (junctions_file == NULL) |
| 1833 |
{ |
| 1834 |
fprintf(stderr, "Error: cannot open BED file %s for writing\n", |
| 1835 |
junctions_file_name.c_str()); |
| 1836 |
exit(1); |
| 1837 |
} |
| 1838 |
|
| 1839 |
FILE* insertions_file = fopen(insertions_file_name.c_str(), "w"); |
| 1840 |
if (insertions_file == NULL) |
| 1841 |
{ |
| 1842 |
fprintf(stderr, "Error: cannot open VCF file %s for writing\n", |
| 1843 |
insertions_file_name.c_str()); |
| 1844 |
exit(1); |
| 1845 |
} |
| 1846 |
|
| 1847 |
FILE* deletions_file = fopen(deletions_file_name.c_str(), "w"); |
| 1848 |
if (deletions_file == NULL) |
| 1849 |
{ |
| 1850 |
fprintf(stderr, "Error: cannot open VCF file %s for writing\n", |
| 1851 |
deletions_file_name.c_str()); |
| 1852 |
exit(1); |
| 1853 |
} |
| 1854 |
|
| 1855 |
FILE* fusions_file = NULL; |
| 1856 |
if (fusion_search) |
| 1857 |
{ |
| 1858 |
fusions_file = fopen(fusions_file_name.c_str(), "w"); |
| 1859 |
if (fusions_file == NULL) |
| 1860 |
{ |
| 1861 |
fprintf(stderr, "Error: cannot open VCF file %s for writing\n", |
| 1862 |
fusions_file_name.c_str()); |
| 1863 |
exit(1); |
| 1864 |
} |
| 1865 |
} |
| 1866 |
|
| 1867 |
driver(accepted_hits_file_name, |
| 1868 |
ref_stream, |
| 1869 |
left_map_filename, |
| 1870 |
left_reads_filename, |
| 1871 |
right_map_filename, |
| 1872 |
right_reads_filename, |
| 1873 |
junctions_file, |
| 1874 |
insertions_file, |
| 1875 |
deletions_file, |
| 1876 |
fusions_file); |
| 1877 |
|
| 1878 |
return 0; |
| 1879 |
} |
