ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/tophat_cpp/fusions.h
Revision: 159
Committed: Mon Jan 30 22:38:08 2012 UTC (7 years, 8 months ago) by gpertea
File size: 5807 byte(s)
Log Message:
wip

Line File contents
1 #ifndef FUSIONS_H
2 #define FUSIONS_H
3 /*
4 * fusions.h
5 * TopHat
6 *
7 * Adapted from junctions.h
8 */
9
10 #include <cstdio>
11 #include <vector>
12 #include <string>
13 #include <set>
14 #include <iostream>
15 #include <fstream>
16 #include <cstring>
17 #include <seqan/sequence.h>
18 #include <seqan/find.h>
19 #include <seqan/file.h>
20
21 #include "bwt_map.h"
22 using namespace std;
23
24 struct Fusion
25 {
26 Fusion (uint32_t ref1, uint32_t ref2, uint32_t l, uint32_t r, uint32_t d = FUSION_FF)
27 : refid1(ref1), refid2(ref2), left(l), right(r), dir(d)
28 {}
29
30 Fusion() : refid1(0), refid2(0), left(0), right(0), dir(FUSION_FF)
31 {}
32
33 uint32_t refid1;
34 uint32_t refid2;
35 uint32_t left;
36 uint32_t right;
37 uint32_t dir;
38
39 bool operator<(const Fusion& rhs) const
40 {
41 if (refid1 < rhs.refid1)
42 return true;
43 else if (refid1 > rhs.refid1)
44 return false;
45
46 if (refid2 < rhs.refid2)
47 return true;
48 else if (refid2 > rhs.refid2)
49 return false;
50
51 if (left < rhs.left)
52 return true;
53 else if (left > rhs.left)
54 return false;
55
56 if (right < rhs.right)
57 return true;
58 else if (right > rhs.right)
59 return false;
60
61 if (dir < rhs.dir)
62 return true;
63 else if(dir > rhs.dir)
64 return false;
65
66 return false;
67 }
68
69 bool operator==(const Fusion& rhs) const
70 {
71 return refid1 == rhs.refid1 &&
72 refid2 == rhs.refid2 &&
73 left == rhs.left &&
74 right == rhs.right &&
75 dir == rhs.dir;
76 }
77 };
78
79 struct FusionPairSupport
80 {
81 FusionPairSupport(int l, int r) :
82 ldist(l),
83 rdist(r)
84 {}
85
86 bool operator<(const FusionPairSupport& rhs) const
87 {
88 if (abs(ldist) + abs(rdist) < abs(rhs.ldist) + abs(rhs.rdist))
89 return true;
90 else
91 return false;
92 }
93
94 int ldist;
95 int rdist;
96 };
97
98 struct FusionSimpleStat
99 {
100 uint32_t count; // # of reads that support the fusion
101 uint32_t edit_dist; // the smallest edit dist among the supporting reads
102 bool skip; //
103 bool left_coincide_with_splice_junction; //
104 bool right_coincide_with_splice_junction; //
105
106 FusionSimpleStat& merge_with(const FusionSimpleStat& other)
107 {
108 if (this == &other)
109 return *this;
110
111 count += other.count;
112 edit_dist = min(edit_dist, other.edit_dist);
113
114 return *this;
115 }
116 };
117
118 struct FusionStat
119 {
120 uint32_t count; // # of reads that support the fusion
121 uint32_t unsupport_count; // # of reads that unsupport the fusion, that is, reads that span one chromosome.
122 uint32_t unsupport_count_pair;
123 uint32_t pair_count; // # of pairs that support the fusion
124 uint32_t pair_count_fusion; // # of pairs that support the fusion where one read spans the fusion
125 float symm;
126
127 string chr1_seq;
128 string chr2_seq;
129
130 static const uint32_t NUM_BASES = 50;
131 vector<uint32_t> left_bases;
132 vector<uint32_t> right_bases;
133
134 vector<uint32_t> diffs;
135
136 uint32_t left_ext;
137 uint32_t right_ext;
138
139 vector<FusionPairSupport> vPairSupport;
140
141 /*
142 * daehwan - this is a metadata indicating whether Fusion is reversed ..
143 * it's true that this is not a good way to ...
144 */
145 bool reversed;
146
147 FusionStat() :
148 left_bases(vector<uint32_t>(NUM_BASES, 0)),
149 right_bases(vector<uint32_t>(NUM_BASES, 0))
150 {
151 count = 0;
152 unsupport_count = 0;
153 unsupport_count_pair = 0;
154 pair_count = 0;
155 pair_count_fusion = 0;
156 symm = 0.0f;
157
158 left_ext = right_ext = 0;
159 reversed = false;
160 }
161
162 FusionStat& merge_with(const FusionStat& other_fusion)
163 {
164 if (this == &other_fusion)
165 return *this;
166
167 count += other_fusion.count;
168 unsupport_count += other_fusion.unsupport_count;
169 unsupport_count_pair += other_fusion.unsupport_count_pair;
170 pair_count += other_fusion.pair_count;
171 pair_count_fusion += other_fusion.pair_count_fusion;
172
173 if (other_fusion.count > 0)
174 {
175 symm = other_fusion.symm;
176 chr1_seq = other_fusion.chr1_seq;
177 chr2_seq = other_fusion.chr2_seq;
178 diffs = other_fusion.diffs;
179 }
180
181 assert (left_bases.size() == right_bases.size());
182 assert (left_bases.size() == other_fusion.left_bases.size());
183 assert (right_bases.size() == other_fusion.right_bases.size());
184 for (size_t i = 0; i < left_bases.size(); ++i)
185 {
186 left_bases[i] += other_fusion.left_bases[i];
187 right_bases[i] += other_fusion.right_bases[i];
188 }
189
190 left_ext = max(left_ext, other_fusion.left_ext);
191 right_ext = max(right_ext, other_fusion.right_ext);
192
193 vPairSupport.insert(vPairSupport.end(), other_fusion.vPairSupport.begin(), other_fusion.vPairSupport.end());
194
195 return *this;
196 }
197 };
198
199 struct fusion_comparison
200 {
201 bool operator()(const Fusion& lhs, const Fusion& rhs)
202 {
203 if (lhs.left != rhs.left)
204 return lhs.left < rhs.left;
205
206 if (lhs.right != rhs.right)
207 return lhs.right < rhs.right;
208
209 if (lhs.refid1 != rhs.refid1)
210 return lhs.refid1 < rhs.refid1;
211
212 if (lhs.refid2 != rhs.refid2)
213 return lhs.refid2 < rhs.refid2;
214
215 if (lhs.dir != rhs.dir)
216 return lhs.dir < rhs.dir;
217
218 return false;
219 }
220 };
221
222 // this is used in segment_juncs
223 typedef std::map<Fusion, FusionSimpleStat> FusionSimpleSet;
224
225 // this is used in tophat_reports
226 typedef std::map<Fusion, FusionStat, fusion_comparison> FusionSet;
227
228 void fusions_from_alignment(const BowtieHit& bh, FusionSet& fusions, RefSequenceTable& rt, bool update_stat = false);
229 void unsupport_fusions(const BowtieHit& bh, FusionSet& fusions, FusionSet& fusions_ref);
230 void print_fusions(FILE* fusions_out, FusionSet& fusions, RefSequenceTable& ref_sequences);
231 void fusions_from_spliced_hit(const BowtieHit& bh, vector<Fusion>& fusions, bool auto_sort = true);
232 void pair_support(const HitsForRead& left_hits, const HitsForRead& right_hits, FusionSet& fusions, FusionSet& fusions_ref);
233
234 void merge_with(FusionSimpleSet& fusions, const FusionSimpleSet& other_fusions);
235 void merge_with(FusionSet& fusions, const FusionSet& other_fusions);
236
237 #endif

Properties

Name Value
svn:executable *