ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/tophat_cpp/junctions.h
Revision: 154
Committed: Tue Jan 24 02:29:21 2012 UTC (7 years, 8 months ago) by gpertea
File size: 3474 byte(s)
Log Message:
massive update with Daehwan's work

Line File contents
1 #ifndef JUNCTIONS_H
2 #define JUNCTIONS_H
3 /*
4 * junctions.h
5 * TopHat
6 *
7 * Created by Cole Trapnell on 11/22/08.
8 * Copyright 2008 Cole Trapnell. All rights reserved.
9 *
10 */
11
12 #include <cstdio>
13 #include <vector>
14 #include <string>
15 #include <set>
16 #include <iostream>
17 #include <fstream>
18 #include <cstring>
19 #include <seqan/sequence.h>
20 #include <seqan/find.h>
21 #include <seqan/file.h>
22
23 #include "bwt_map.h"
24
25 using namespace std;
26
27 struct Junction
28 {
29 Junction (uint32_t ref, uint32_t l, uint32_t r, bool a, int sc = 0)
30 : refid(ref), left(l), right(r), antisense(a), skip_count(sc){}
31 Junction() : refid(0), left(0), right(0), antisense(false), skip_count(0) {}
32 uint32_t refid;
33 uint32_t left;
34 uint32_t right;
35 bool antisense;
36
37 int skip_count;
38
39 bool operator<(const Junction& rhs) const
40 {
41 if (refid < rhs.refid)
42 return true;
43 else if (refid > rhs.refid)
44 return false;
45
46 if (left < rhs.left)
47 return true;
48 else if (left > rhs.left)
49 return false;
50
51 if (right < rhs.right)
52 return true;
53 else if (right > rhs.right)
54 return false;
55
56 return antisense < rhs.antisense;
57 }
58
59 bool operator==(const Junction& rhs) const
60 {
61 return (refid == rhs.refid && left == rhs.left && right == rhs.right && antisense == rhs.antisense);
62 }
63
64 #if !NDEBUG
65 bool valid() const
66 {
67 return refid != VMAXINT32 && left < right && (left != right);
68 }
69 #endif
70 };
71
72 struct skip_count_lt
73 {
74 bool operator()(const Junction& lhs, const Junction& rhs)
75 {
76 if (lhs.skip_count != rhs.skip_count)
77 return lhs.skip_count < rhs.skip_count;
78 return lhs < rhs;
79 }
80 };
81
82 struct JunctionStats
83 {
84 JunctionStats() : left_extent(0), right_extent(0), left_exon_doc(0), right_exon_doc(0),
85 min_splice_mms(0), supporting_hits(0), gtf_match(false), accepted(false) {}
86
87 JunctionStats& merge_with(const JunctionStats& other)
88 {
89 if (this == &other)
90 return *this;
91
92 left_extent = max(left_extent, other.left_extent);
93 right_extent = max(right_extent, other.right_extent);
94 min_splice_mms = min(min_splice_mms, other.min_splice_mms);
95 supporting_hits += other.supporting_hits;
96
97 return *this;
98 }
99
100 int left_extent;
101 int right_extent;
102 int left_exon_doc;
103 int right_exon_doc;
104 int min_splice_mms;
105 int supporting_hits;
106 bool gtf_match;
107 bool accepted;
108 };
109
110 typedef std::map<Junction, JunctionStats> JunctionSet;
111
112 // This routine DOES NOT set the real refid!
113 pair<Junction, JunctionStats> junction_from_spliced_hit(const BowtieHit& h);
114
115 void print_junction(FILE* junctions_out,
116 const string& name,
117 const Junction& j,
118 const JunctionStats& s,
119 uint32_t junc_id);
120
121
122 void junctions_from_alignment(const BowtieHit& spliced_alignment,
123 JunctionSet& junctions);
124
125 void accept_valid_junctions(JunctionSet& junctions,
126 const uint32_t refid,
127 const vector<unsigned short>& DoC,
128 double min_isoform_fraction);
129
130 void accept_all_junctions(JunctionSet& junctions,
131 const uint32_t refid);
132
133 void print_junctions(FILE* junctions_out,
134 const JunctionSet& junctions,
135 RefSequenceTable& ref_sequences);
136
137 bool accept_if_valid(const Junction& j, JunctionStats& s);
138
139 void filter_junctions(JunctionSet& junctions, const JunctionSet& gtf_junctions);
140
141 void get_junctions_from_hits(HitStream& hit_stream,
142 ReadTable& it,
143 JunctionSet& junctions);
144
145 void merge_with(JunctionSet& juncs, const JunctionSet& other_juncs);
146
147 #endif