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 (8 years, 8 months ago) by gpertea
File size: 5807 byte(s)
Log Message:
wip

Line User Rev File contents
1 gpertea 158 #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 gpertea 159 symm = other_fusion.symm;
176 gpertea 158 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 *