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"
#include "junctions.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),
    fusion(true),
    inner_dist(99999999),
    alignment_score(std::numeric_limits<int>::min())
  {}

InsertAlignmentGrade(const BowtieHit& h1, bool fusion = false) : 
  too_close(false),
    too_far(false),
    num_spliced(0),
    num_mapped(0),
    opposite_strands(false),
    consistent_splices(false),
    edit_dist(0x1F),
    fusion(fusion),
    num_alignments(0),
    inner_dist(99999999),
    alignment_score(std::numeric_limits<int>::min())
  {
    if (!h1.contiguous())
      num_spliced++;
    
    num_mapped = 1;
    
    longest_ref_skip = min(0x3FFFFu, (unsigned int)get_longest_ref_skip(h1) / 100);
    edit_dist = h1.edit_dist();
    num_alignments = 1;
  }
  
InsertAlignmentGrade(const BowtieHit& h1, 
                     const BowtieHit& h2,
                     const JunctionSet& junctions,
                     bool fusion = false) :
  too_close(false),
    too_far(false),
    num_spliced(0),
    num_mapped(0),
    opposite_strands(false),
    consistent_splices(false),
    edit_dist(0x1F),
    fusion(fusion),
    num_alignments(0),
    alignment_score(std::numeric_limits<int>::min())
  {
    int min_inner_distance = inner_dist_mean - inner_dist_std_dev;
    int max_inner_distance = inner_dist_mean + inner_dist_std_dev;
  
    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));
    longest_ref_skip = min (ls / 100, 0x3FFFFu);
    edit_dist = h1.edit_dist() + h2.edit_dist();
    num_alignments = 1;
    assert(!(too_far && too_close));

    if (too_far && !fusion)
      {
        int inner1, inner2;
        if (h1.left() < h2.left())
          inner1 = h1.right(), inner2 = h2.left();
        else
          inner1 = h2.right(), inner2 = h1.left();
        
        JunctionSet::const_iterator lb, ub;

        lb = junctions.upper_bound(Junction(h1.ref_id(), inner1, inner1, true));
        ub = junctions.lower_bound(Junction(h1.ref_id(), inner2, inner2, false));
        while (lb != ub && lb != junctions.end())
          {
            const Junction& junction = lb->first;
            const JunctionStats& junction_stat = lb->second;
            if (inner1 <= (int)junction.left && inner2 >= (int)junction.right &&
                junction_stat.supporting_hits >= 10)
              {
                int temp_dist = junction.left - inner1 + inner2 - junction.right;
                if (temp_dist >= min_inner_distance && temp_dist <= max_inner_distance)
                  {
                    // daehwan - for debugging purposes
                    /*
                    fprintf(stderr, "read id: %d\t %s %d-%d %d-%d %s\t %d->%d) \n",
                            h1.insert_id(), print_cigar(h1.cigar()).c_str(), h1.left(), h1.right(),
                            h2.left(), h2.right(), print_cigar(h2.cigar()).c_str(),
                            inner_dist, temp_dist);
                    */
                    
                    inner_dist = temp_dist;
                    too_far = false;

                    break;
                  }
              }
            
            ++lb;
          }
      }

    //
    static const int penalty_for_long_inner_dist = bowtie2_max_penalty;
    alignment_score = h1.alignment_score() + h2.alignment_score();
    if (!fusion)
      {
        if (too_far)
          {
            int penalty = penalty_for_long_inner_dist;
            if (inner_dist - max_inner_distance < inner_dist_std_dev)
              {
                penalty = penalty_for_long_inner_dist / 2;
              }
              
            alignment_score -= penalty;
          }
        else if (too_close)
          {
            int penalty = min(penalty_for_long_inner_dist/2,
                              (min_inner_distance - inner_dist) / inner_dist_std_dev + 1);
            alignment_score -= penalty;
          }
        
        static const int penalty_for_same_strand = bowtie2_max_penalty;
        if (!opposite_strands)
          alignment_score -= penalty_for_same_strand;
      }
  }
  
  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;
    fusion = rhs.fusion;
    num_alignments = rhs.num_alignments;
    inner_dist = rhs.inner_dist;
    alignment_score = rhs.alignment_score;
    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_element(gaps.begin(), gaps.end(), gap_lt);
    return abs(max_itr->second - max_itr->first);
  }
  
  // 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;
  }

  int align_score() { return alignment_score; }

  bool is_fusion() const { return fusion; }
  
  bool too_close;
  bool too_far;
  
  uint8_t num_spliced;
  uint8_t num_mapped;
  
  bool opposite_strands;
  bool consistent_splices;
  uint32_t longest_ref_skip; // in 100s of bp
  unsigned char edit_dist;
  bool fusion;
  int num_alignments; // number of equally good alignments for the insert 
  int inner_dist; // distance between inner edges of mates

  int alignment_score;
};              

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:	245
Committed:	Wed May 16 21:51:01 2012 UTC (12 years, 3 months ago) by gpertea
File size:	7167 byte(s)
Log Message:	merge-sync attempt
Line	File contents
1	#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	#include "junctions.h"
14
15	struct InsertAlignment
16	{
17	InsertAlignment(uint64_t _refid,
18	BowtieHit* _left_alignment,
19	BowtieHit* _right_alignment) :
20	refid(_refid),
21	left_alignment(_left_alignment),
22	right_alignment(_right_alignment) {}
23
24	uint64_t refid;
25	BowtieHit* left_alignment;
26	BowtieHit* right_alignment;
27	};
28
29	pair<int, int> pair_distances(const BowtieHit& h1, const BowtieHit& h2);
30
31	bool gap_lt(const pair<int, int>& lhs, const pair<int, int>& rhs);
32
33	struct InsertAlignmentGrade
34	{
35	InsertAlignmentGrade() :
36	too_close(false),
37	too_far(false),
38	num_spliced(0),
39	num_mapped(0),
40	opposite_strands(false),
41	consistent_splices(false),
42	longest_ref_skip(0x7FFFFu),
43	edit_dist(0x1F),
44	fusion(true),
45	inner_dist(99999999),
46	alignment_score(std::numeric_limits<int>::min())
47	{}
48
49	InsertAlignmentGrade(const BowtieHit& h1, bool fusion = false) :
50	too_close(false),
51	too_far(false),
52	num_spliced(0),
53	num_mapped(0),
54	opposite_strands(false),
55	consistent_splices(false),
56	edit_dist(0x1F),
57	fusion(fusion),
58	num_alignments(0),
59	inner_dist(99999999),
60	alignment_score(std::numeric_limits<int>::min())
61	{
62	if (!h1.contiguous())
63	num_spliced++;
64
65	num_mapped = 1;
66
67	longest_ref_skip = min(0x3FFFFu, (unsigned int)get_longest_ref_skip(h1) / 100);
68	edit_dist = h1.edit_dist();
69	num_alignments = 1;
70	}
71
72	InsertAlignmentGrade(const BowtieHit& h1,
73	const BowtieHit& h2,
74	const JunctionSet& junctions,
75	bool fusion = false) :
76	too_close(false),
77	too_far(false),
78	num_spliced(0),
79	num_mapped(0),
80	opposite_strands(false),
81	consistent_splices(false),
82	edit_dist(0x1F),
83	fusion(fusion),
84	num_alignments(0),
85	alignment_score(std::numeric_limits<int>::min())
86	{
87	int min_inner_distance = inner_dist_mean - inner_dist_std_dev;
88	int max_inner_distance = inner_dist_mean + inner_dist_std_dev;
89
90	pair<int, int> distances = pair_distances(h1, h2);
91	inner_dist = distances.second;
92
93	num_mapped = 2;
94
95	if (!h1.contiguous())
96	num_spliced++;
97	if (!h2.contiguous())
98	num_spliced++;
99
100	too_far = (inner_dist > max_inner_distance);
101	too_close = (inner_dist < min_inner_distance);
102	opposite_strands = (h1.antisense_align() != h2.antisense_align());
103	consistent_splices = (num_spliced == 2 && h1.antisense_splice() == h2.antisense_splice());
104
105	uint32_t ls = max(get_longest_ref_skip(h1), get_longest_ref_skip(h2));
106	longest_ref_skip = min (ls / 100, 0x3FFFFu);
107	edit_dist = h1.edit_dist() + h2.edit_dist();
108	num_alignments = 1;
109	assert(!(too_far && too_close));
110
111	if (too_far && !fusion)
112	{
113	int inner1, inner2;
114	if (h1.left() < h2.left())
115	inner1 = h1.right(), inner2 = h2.left();
116	else
117	inner1 = h2.right(), inner2 = h1.left();
118
119	JunctionSet::const_iterator lb, ub;
120
121	lb = junctions.upper_bound(Junction(h1.ref_id(), inner1, inner1, true));
122	ub = junctions.lower_bound(Junction(h1.ref_id(), inner2, inner2, false));
123	while (lb != ub && lb != junctions.end())
124	{
125	const Junction& junction = lb->first;
126	const JunctionStats& junction_stat = lb->second;
127	if (inner1 <= (int)junction.left && inner2 >= (int)junction.right &&
128	junction_stat.supporting_hits >= 10)
129	{
130	int temp_dist = junction.left - inner1 + inner2 - junction.right;
131	if (temp_dist >= min_inner_distance && temp_dist <= max_inner_distance)
132	{
133	// daehwan - for debugging purposes
134	/*
135	fprintf(stderr, "read id: %d\t %s %d-%d %d-%d %s\t %d->%d) \n",
136	h1.insert_id(), print_cigar(h1.cigar()).c_str(), h1.left(), h1.right(),
137	h2.left(), h2.right(), print_cigar(h2.cigar()).c_str(),
138	inner_dist, temp_dist);
139	*/
140
141	inner_dist = temp_dist;
142	too_far = false;
143
144	break;
145	}
146	}
147
148	++lb;
149	}
150	}
151
152	//
153	static const int penalty_for_long_inner_dist = bowtie2_max_penalty;
154	alignment_score = h1.alignment_score() + h2.alignment_score();
155	if (!fusion)
156	{
157	if (too_far)
158	{
159	int penalty = penalty_for_long_inner_dist;
160	if (inner_dist - max_inner_distance < inner_dist_std_dev)
161	{
162	penalty = penalty_for_long_inner_dist / 2;
163	}
164
165	alignment_score -= penalty;
166	}
167	else if (too_close)
168	{
169	int penalty = min(penalty_for_long_inner_dist/2,
170	(min_inner_distance - inner_dist) / inner_dist_std_dev + 1);
171	alignment_score -= penalty;
172	}
173
174	static const int penalty_for_same_strand = bowtie2_max_penalty;
175	if (!opposite_strands)
176	alignment_score -= penalty_for_same_strand;
177	}
178	}
179
180	InsertAlignmentGrade& operator=(const InsertAlignmentGrade& rhs)
181	{
182	too_close = rhs.too_close;
183	too_far = rhs.too_far;
184	num_spliced = rhs.num_spliced;
185	num_mapped = rhs.num_mapped;
186	opposite_strands = rhs.opposite_strands;
187	consistent_splices = rhs.consistent_splices;
188	longest_ref_skip = rhs.longest_ref_skip;
189	edit_dist = rhs.edit_dist;
190	fusion = rhs.fusion;
191	num_alignments = rhs.num_alignments;
192	inner_dist = rhs.inner_dist;
193	alignment_score = rhs.alignment_score;
194	return *this;
195	}
196
197	static int get_longest_ref_skip(const BowtieHit& h1)
198	{
199	vector<pair<int, int> > gaps;
200	h1.gaps(gaps);
201	if (gaps.empty())
202	return 0;
203
204	vector<pair<int, int> >::iterator max_itr = max_element(gaps.begin(), gaps.end(), gap_lt);
205	return abs(max_itr->second - max_itr->first);
206	}
207
208	// Returns true if rhs is a "happier" alignment for the ends of this insert
209	// than this InsertStatus.
210	bool operator<(const InsertAlignmentGrade& rhs);
211
212	bool happy() const
213	{
214	return num_mapped == 2 && opposite_strands && (num_spliced != 2 \|\| consistent_splices) && !too_far;
215	}
216
217	int align_score() { return alignment_score; }
218
219	bool is_fusion() const { return fusion; }
220
221	bool too_close;
222	bool too_far;
223
224	uint8_t num_spliced;
225	uint8_t num_mapped;
226
227	bool opposite_strands;
228	bool consistent_splices;
229	uint32_t longest_ref_skip; // in 100s of bp
230	unsigned char edit_dist;
231	bool fusion;
232	int num_alignments; // number of equally good alignments for the insert
233	int inner_dist; // distance between inner edges of mates
234
235	int alignment_score;
236	};
237
238	typedef vector<pair<InsertAlignmentGrade, vector<InsertAlignment> > > BestInsertAlignmentTable;
239	void accept_valid_hits(BestInsertAlignmentTable& best_status_for_inserts);
240	void accept_all_best_hits(BestInsertAlignmentTable& best_status_for_inserts);
241
242	void best_insert_mappings(uint64_t refid,
243	ReadTable& it,
244	HitList& hits1_in_ref,
245	HitList& hits2_in_ref,
246	BestInsertAlignmentTable& best_insert_alignments,
247	bool prefer_shorter_pairs = false);
248
249
250	void insert_best_pairings(RefSequenceTable& rt,
251	ReadTable& it,
252	HitTable& hits1,
253	HitTable& hits2,
254	BestInsertAlignmentTable& best_pairings,
255	bool prefer_shorter_pairs = false);
256	#endif