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

Line File contents
1 /*
2 * align_status.cpp
3 * TopHat
4 *
5 * Created by Ryan Kelley on 11/09/2010.
6 *
7 */
8
9 #ifdef HAVE_CONFIG_H
10 #include <config.h>
11 #endif
12
13 #include <cassert>
14 #include <cstdlib>
15 #include <iostream>
16 #include <set>
17 #include <vector>
18
19 #include "common.h"
20 #include "bwt_map.h"
21 #include "tokenize.h"
22 #include "reads.h"
23 #include "align_status.h"
24
25 using namespace std;
26
27 AlignStatus::AlignStatus()
28 {
29 _aligned = false;
30 _indelFreeAlignment = false;
31 _unannotatedSpliceFreeAlignment = false;
32 _edit_dist = 0xFF;
33 _fusionFreeAlignment = 0;
34 }
35
36 /**
37 * Parse the cigar string of a BowtieHit in order to determine the alignment status.
38 */
39 AlignStatus::AlignStatus(const BowtieHit& bh, const JunctionSet& gtf_junctions){
40 const vector<CigarOp>& cigar = bh.cigar();
41 _aligned = cigar.size() > 0;
42 _indelFreeAlignment = true;
43 _fusionFreeAlignment = true;
44 _unannotatedSpliceFreeAlignment = true;
45 _edit_dist = bh.edit_dist();
46 int j = bh.left();
47 for (size_t c = 0 ; c < cigar.size(); ++c)
48 {
49 Junction junc;
50 switch(cigar[c].opcode)
51 {
52 case REF_SKIP:
53 junc.refid = bh.ref_id();
54 junc.left = j;
55 junc.right = junc.left + cigar[c].length;
56 junc.antisense = bh.antisense_splice();
57 j += cigar[c].length;
58
59 if (gtf_junctions.find(junc) == gtf_junctions.end())
60 _unannotatedSpliceFreeAlignment = false;
61 break;
62 case MATCH:
63 j += cigar[c].length;
64 break;
65 case DEL:
66 j += cigar[c].length;
67 _indelFreeAlignment = false;
68 break;
69 case INS:
70 _indelFreeAlignment = false;
71 break;
72 case FUSION_FF:
73 case FUSION_FR:
74 case FUSION_RF:
75 case FUSION_RR:
76 _fusionFreeAlignment = false;
77 break;
78 default:
79 break;
80 }
81 }
82 }
83
84 /**
85 * Establish an ordering on alignments.
86 * Prefer aligned reads over unaligned reads
87 * Within the space of aligned reads
88 * prefer splice-free reads over splice reads, and
89 * indel-free reads over indel reads.
90 * If a read can either be indel-free or splice-free,
91 * prefer the indel-free alignment
92 */
93 bool AlignStatus::operator<(const AlignStatus& rhs) const
94 {
95 if (rhs._aligned != _aligned) return rhs._aligned;
96 if (rhs._edit_dist!=_edit_dist)
97 return rhs._edit_dist < _edit_dist;
98
99 // int lhs_value = _aligned ? 1 : 0;
100 int lhs_value = _fusionFreeAlignment ? 4 : 0;
101 lhs_value += _indelFreeAlignment ? 4 : 0;
102 lhs_value += _unannotatedSpliceFreeAlignment ? 2 : 0;
103
104 // int rhs_value = rhs._aligned ? 1 : 0;
105 int rhs_value = rhs._fusionFreeAlignment ? 4 : 0;
106 rhs_value += rhs._indelFreeAlignment ? 4 : 0;
107 rhs_value += rhs._unannotatedSpliceFreeAlignment ? 2 : 0;
108
109 return lhs_value < rhs_value;
110 }
111
112 /**
113 * Alignments are only equal if all fields are identical.
114 */
115 bool AlignStatus::operator==(const AlignStatus& rhs) const
116 {
117 return ((_aligned == rhs._aligned) && (rhs._edit_dist ==_edit_dist) &&
118 (_indelFreeAlignment == rhs._indelFreeAlignment) &&
119 (_unannotatedSpliceFreeAlignment == rhs._unannotatedSpliceFreeAlignment) &&
120 (_fusionFreeAlignment == rhs._fusionFreeAlignment));
121 }
122
123 bool AlignStatus::operator!=(const AlignStatus& rhs) const
124 {
125 return !((*this) == rhs);
126 }