gclib/tophat_cpp/inserts.h

#ifndef INSERTS_H
#define INSERTS_H
/*
 *  inserts.h
 *  TopHat
 *
 *  Created by Cole Trapnell on 1/14/09.
 *  Copyright 2009 Cole Trapnell. All rights reserved.
 *
 */

#include "bwt_map.h"

struct InsertAlignment
{
        InsertAlignment(uint64_t _refid, 
                                        BowtieHit* _left_alignment, 
                                        BowtieHit* _right_alignment) : 
        refid(_refid), 
        left_alignment(_left_alignment),
        right_alignment(_right_alignment) {}
        
        uint64_t refid;
        BowtieHit* left_alignment;
        BowtieHit* right_alignment;
};

pair<int, int> pair_distances(const BowtieHit& h1, const BowtieHit& h2);

bool gap_lt(const pair<int, int>& lhs, const pair<int, int>& rhs);

struct InsertAlignmentGrade
{
        InsertAlignmentGrade() : 
                too_close(false),
                too_far(false),
                num_spliced(0),
                num_mapped(0),
                opposite_strands(false),
                consistent_splices(false),
                longest_ref_skip(0x7FFFFu),
                edit_dist(0x1F), 
                inner_dist(99999999){}
        
        InsertAlignmentGrade(const BowtieHit& h1) : 
                too_close(false),
                too_far(false),
                num_spliced(0),
                num_mapped(0),
                opposite_strands(false),
                consistent_splices(false),
                edit_dist(0x1F),
                num_alignments(0),
                inner_dist(99999999)
        {
                if (!h1.contiguous())
                        num_spliced++;
                num_mapped = 1;
                
                longest_ref_skip = min(0x7FFFFu, (unsigned int)get_longest_ref_skip(h1));
                edit_dist = h1.edit_dist();
                
                longest_ref_skip /= 100;
                num_alignments = 1;
                
        }
        
        InsertAlignmentGrade(const BowtieHit& h1, 
                             const BowtieHit& h2, 
                             int min_inner_distance,
                             int max_inner_distance) :
                  too_close(false),
                  too_far(false),
                  num_spliced(0),
                  num_mapped(0),
                  opposite_strands(false),
                  consistent_splices(false),
                  edit_dist(0x1F),
                  num_alignments(0)
        {
                pair<int, int> distances = pair_distances(h1,h2);
                inner_dist = distances.second;
                
                num_mapped = 2;
                
                if(!h1.contiguous())
                        num_spliced++;
                if(!h2.contiguous())
                        num_spliced++;
                
                too_far = (inner_dist > max_inner_distance);
                
                too_close = (inner_dist < min_inner_distance);
                
                opposite_strands = (h1.antisense_align() != h2.antisense_align());
                
                consistent_splices = (num_spliced == 2 &&
                                                          h1.antisense_splice() == h2.antisense_splice());
                
                uint32_t ls = max(get_longest_ref_skip(h1), get_longest_ref_skip(h2));
                ls /= 100;
                
                longest_ref_skip = min (ls, 0x7FFFFu);
                
                edit_dist = h1.edit_dist() + h2.edit_dist();
                
                num_alignments = 1;
                
                assert(!(too_far && too_close));
        }
        
        InsertAlignmentGrade& operator=(const InsertAlignmentGrade& rhs)
        {
                too_close = rhs.too_close;
                too_far = rhs.too_far;
                num_spliced = rhs.num_spliced;
                num_mapped = rhs.num_mapped;
                opposite_strands = rhs.opposite_strands;
                consistent_splices = rhs.consistent_splices;
                longest_ref_skip = rhs.longest_ref_skip;
                edit_dist = rhs.edit_dist;
                num_alignments = rhs.num_alignments;
                inner_dist = rhs.inner_dist;
                return *this;
        }
        
        static int get_longest_ref_skip(const BowtieHit& h1)
        {
                vector<pair<int, int> > gaps;
                h1.gaps(gaps);
                if (gaps.empty())
                {
                        return 0;
                }
                vector<pair<int, int> >::iterator max_itr;
                max_itr = max_element(gaps.begin(), gaps.end(), gap_lt);
                return min(max_itr->second - max_itr->first, 0x7FFFF);
        }
        
        // Returns true if rhs is a "happier" alignment for the ends of this insert
        // than this InsertStatus.
        bool operator<(const InsertAlignmentGrade& rhs);
        
        bool happy() const
        {
                return num_mapped == 2 && opposite_strands && (num_spliced != 2 || consistent_splices) && !too_far;
        }
        
        bool too_close : 1;
        bool too_far : 1;
        
        uint8_t num_spliced : 2;
        
        uint8_t num_mapped : 2;
        
        bool opposite_strands : 1;
        bool consistent_splices : 1;
        uint32_t longest_ref_skip : 19; // in 100s of bp
        unsigned char edit_dist : 5;
        int num_alignments; // number of equally good alignments for the insert 
        int inner_dist; // distance between inner edges of mates
};              

typedef vector<pair<InsertAlignmentGrade, vector<InsertAlignment> > > BestInsertAlignmentTable;
void accept_valid_hits(BestInsertAlignmentTable& best_status_for_inserts);
void accept_all_best_hits(BestInsertAlignmentTable& best_status_for_inserts);

void best_insert_mappings(uint64_t refid,
                                                  ReadTable& it,
                                                  HitList& hits1_in_ref,
                                                  HitList& hits2_in_ref,
                                                  BestInsertAlignmentTable& best_insert_alignments,
                                                  bool prefer_shorter_pairs = false);


void insert_best_pairings(RefSequenceTable& rt,
                                                  ReadTable& it,
                          HitTable& hits1,
                          HitTable& hits2,
                          BestInsertAlignmentTable& best_pairings, 
                                                  bool prefer_shorter_pairs = false);
#endif
Revision:	29
Committed:	Tue Aug 2 21:24:54 2011 UTC (13 years, 1 month ago) by gpertea
File size:	4590 byte(s)
Log Message:	adding tophat source work
Line	User	Rev	File contents
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