1 |
gpertea |
29 |
#ifndef INSERTS_H |
2 |
|
|
#define INSERTS_H |
3 |
|
|
/* |
4 |
|
|
* inserts.h |
5 |
|
|
* TopHat |
6 |
|
|
* |
7 |
|
|
* Created by Cole Trapnell on 1/14/09. |
8 |
|
|
* Copyright 2009 Cole Trapnell. All rights reserved. |
9 |
|
|
* |
10 |
|
|
*/ |
11 |
|
|
|
12 |
|
|
#include "bwt_map.h" |
13 |
|
|
|
14 |
|
|
struct InsertAlignment |
15 |
|
|
{ |
16 |
|
|
InsertAlignment(uint64_t _refid, |
17 |
|
|
BowtieHit* _left_alignment, |
18 |
|
|
BowtieHit* _right_alignment) : |
19 |
|
|
refid(_refid), |
20 |
|
|
left_alignment(_left_alignment), |
21 |
|
|
right_alignment(_right_alignment) {} |
22 |
|
|
|
23 |
|
|
uint64_t refid; |
24 |
|
|
BowtieHit* left_alignment; |
25 |
|
|
BowtieHit* right_alignment; |
26 |
|
|
}; |
27 |
|
|
|
28 |
|
|
pair<int, int> pair_distances(const BowtieHit& h1, const BowtieHit& h2); |
29 |
|
|
|
30 |
|
|
bool gap_lt(const pair<int, int>& lhs, const pair<int, int>& rhs); |
31 |
|
|
|
32 |
|
|
struct InsertAlignmentGrade |
33 |
|
|
{ |
34 |
|
|
InsertAlignmentGrade() : |
35 |
|
|
too_close(false), |
36 |
|
|
too_far(false), |
37 |
|
|
num_spliced(0), |
38 |
|
|
num_mapped(0), |
39 |
|
|
opposite_strands(false), |
40 |
|
|
consistent_splices(false), |
41 |
|
|
longest_ref_skip(0x7FFFFu), |
42 |
|
|
edit_dist(0x1F), |
43 |
|
|
inner_dist(99999999){} |
44 |
|
|
|
45 |
|
|
InsertAlignmentGrade(const BowtieHit& h1) : |
46 |
|
|
too_close(false), |
47 |
|
|
too_far(false), |
48 |
|
|
num_spliced(0), |
49 |
|
|
num_mapped(0), |
50 |
|
|
opposite_strands(false), |
51 |
|
|
consistent_splices(false), |
52 |
|
|
edit_dist(0x1F), |
53 |
|
|
num_alignments(0), |
54 |
|
|
inner_dist(99999999) |
55 |
|
|
{ |
56 |
|
|
if (!h1.contiguous()) |
57 |
|
|
num_spliced++; |
58 |
|
|
num_mapped = 1; |
59 |
|
|
|
60 |
|
|
longest_ref_skip = min(0x7FFFFu, (unsigned int)get_longest_ref_skip(h1)); |
61 |
|
|
edit_dist = h1.edit_dist(); |
62 |
|
|
|
63 |
|
|
longest_ref_skip /= 100; |
64 |
|
|
num_alignments = 1; |
65 |
|
|
|
66 |
|
|
} |
67 |
|
|
|
68 |
|
|
InsertAlignmentGrade(const BowtieHit& h1, |
69 |
|
|
const BowtieHit& h2, |
70 |
|
|
int min_inner_distance, |
71 |
|
|
int max_inner_distance) : |
72 |
|
|
too_close(false), |
73 |
|
|
too_far(false), |
74 |
|
|
num_spliced(0), |
75 |
|
|
num_mapped(0), |
76 |
|
|
opposite_strands(false), |
77 |
|
|
consistent_splices(false), |
78 |
|
|
edit_dist(0x1F), |
79 |
|
|
num_alignments(0) |
80 |
|
|
{ |
81 |
|
|
pair<int, int> distances = pair_distances(h1,h2); |
82 |
|
|
inner_dist = distances.second; |
83 |
|
|
|
84 |
|
|
num_mapped = 2; |
85 |
|
|
|
86 |
|
|
if(!h1.contiguous()) |
87 |
|
|
num_spliced++; |
88 |
|
|
if(!h2.contiguous()) |
89 |
|
|
num_spliced++; |
90 |
|
|
|
91 |
|
|
too_far = (inner_dist > max_inner_distance); |
92 |
|
|
|
93 |
|
|
too_close = (inner_dist < min_inner_distance); |
94 |
|
|
|
95 |
|
|
opposite_strands = (h1.antisense_align() != h2.antisense_align()); |
96 |
|
|
|
97 |
|
|
consistent_splices = (num_spliced == 2 && |
98 |
|
|
h1.antisense_splice() == h2.antisense_splice()); |
99 |
|
|
|
100 |
|
|
uint32_t ls = max(get_longest_ref_skip(h1), get_longest_ref_skip(h2)); |
101 |
|
|
ls /= 100; |
102 |
|
|
|
103 |
|
|
longest_ref_skip = min (ls, 0x7FFFFu); |
104 |
|
|
|
105 |
|
|
edit_dist = h1.edit_dist() + h2.edit_dist(); |
106 |
|
|
|
107 |
|
|
num_alignments = 1; |
108 |
|
|
|
109 |
|
|
assert(!(too_far && too_close)); |
110 |
|
|
} |
111 |
|
|
|
112 |
|
|
InsertAlignmentGrade& operator=(const InsertAlignmentGrade& rhs) |
113 |
|
|
{ |
114 |
|
|
too_close = rhs.too_close; |
115 |
|
|
too_far = rhs.too_far; |
116 |
|
|
num_spliced = rhs.num_spliced; |
117 |
|
|
num_mapped = rhs.num_mapped; |
118 |
|
|
opposite_strands = rhs.opposite_strands; |
119 |
|
|
consistent_splices = rhs.consistent_splices; |
120 |
|
|
longest_ref_skip = rhs.longest_ref_skip; |
121 |
|
|
edit_dist = rhs.edit_dist; |
122 |
|
|
num_alignments = rhs.num_alignments; |
123 |
|
|
inner_dist = rhs.inner_dist; |
124 |
|
|
return *this; |
125 |
|
|
} |
126 |
|
|
|
127 |
|
|
static int get_longest_ref_skip(const BowtieHit& h1) |
128 |
|
|
{ |
129 |
|
|
vector<pair<int, int> > gaps; |
130 |
|
|
h1.gaps(gaps); |
131 |
|
|
if (gaps.empty()) |
132 |
|
|
{ |
133 |
|
|
return 0; |
134 |
|
|
} |
135 |
|
|
vector<pair<int, int> >::iterator max_itr; |
136 |
|
|
max_itr = max_element(gaps.begin(), gaps.end(), gap_lt); |
137 |
|
|
return min(max_itr->second - max_itr->first, 0x7FFFF); |
138 |
|
|
} |
139 |
|
|
|
140 |
|
|
// Returns true if rhs is a "happier" alignment for the ends of this insert |
141 |
|
|
// than this InsertStatus. |
142 |
|
|
bool operator<(const InsertAlignmentGrade& rhs); |
143 |
|
|
|
144 |
|
|
bool happy() const |
145 |
|
|
{ |
146 |
|
|
return num_mapped == 2 && opposite_strands && (num_spliced != 2 || consistent_splices) && !too_far; |
147 |
|
|
} |
148 |
|
|
|
149 |
|
|
bool too_close : 1; |
150 |
|
|
bool too_far : 1; |
151 |
|
|
|
152 |
|
|
uint8_t num_spliced : 2; |
153 |
|
|
|
154 |
|
|
uint8_t num_mapped : 2; |
155 |
|
|
|
156 |
|
|
bool opposite_strands : 1; |
157 |
|
|
bool consistent_splices : 1; |
158 |
|
|
uint32_t longest_ref_skip : 19; // in 100s of bp |
159 |
|
|
unsigned char edit_dist : 5; |
160 |
|
|
int num_alignments; // number of equally good alignments for the insert |
161 |
|
|
int inner_dist; // distance between inner edges of mates |
162 |
|
|
}; |
163 |
|
|
|
164 |
|
|
typedef vector<pair<InsertAlignmentGrade, vector<InsertAlignment> > > BestInsertAlignmentTable; |
165 |
|
|
void accept_valid_hits(BestInsertAlignmentTable& best_status_for_inserts); |
166 |
|
|
void accept_all_best_hits(BestInsertAlignmentTable& best_status_for_inserts); |
167 |
|
|
|
168 |
|
|
void best_insert_mappings(uint64_t refid, |
169 |
|
|
ReadTable& it, |
170 |
|
|
HitList& hits1_in_ref, |
171 |
|
|
HitList& hits2_in_ref, |
172 |
|
|
BestInsertAlignmentTable& best_insert_alignments, |
173 |
|
|
bool prefer_shorter_pairs = false); |
174 |
|
|
|
175 |
|
|
|
176 |
|
|
void insert_best_pairings(RefSequenceTable& rt, |
177 |
|
|
ReadTable& it, |
178 |
|
|
HitTable& hits1, |
179 |
|
|
HitTable& hits2, |
180 |
|
|
BestInsertAlignmentTable& best_pairings, |
181 |
|
|
bool prefer_shorter_pairs = false); |
182 |
|
|
#endif |