ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/tophat_cpp/align_status.cpp
Revision: 135
Committed: Mon Dec 12 22:28:38 2011 UTC (8 years, 5 months ago) by gpertea
File size: 2738 byte(s)
Log Message:
wip - SplicedSAMHitFactory() still not implemented

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 }
34
35 /**
36 * Parse the cigar string of a BowtieHit in order to determine the alignment status.
37 */
38 AlignStatus::AlignStatus(const BowtieHit& bh, const JunctionSet& gtf_junctions)
39 {
40 const vector<CigarOp>& cigar = bh.cigar();
41 _aligned = cigar.size() > 0;
42 _indelFreeAlignment = true;
43 _unannotatedSpliceFreeAlignment = true;
44 _edit_dist = bh.edit_dist();
45 int j = bh.left();
46 for (size_t c = 0 ; c < cigar.size(); ++c)
47 {
48 Junction junc;
49 switch(cigar[c].opcode)
50 {
51 case REF_SKIP:
52 junc.refid = bh.ref_id();
53 junc.left = j;
54 junc.right = junc.left + cigar[c].length;
55 junc.antisense = bh.antisense_splice();
56 j += cigar[c].length;
57
58 if (gtf_junctions.find(junc) == gtf_junctions.end())
59 _unannotatedSpliceFreeAlignment = false;
60 break;
61 case MATCH:
62 j += cigar[c].length;
63 break;
64 case DEL:
65 j += cigar[c].length;
66 _indelFreeAlignment = false;
67 break;
68 case INS:
69 _indelFreeAlignment = false;
70 break;
71 default:
72 break;
73 }
74 }
75 }
76
77 /**
78 * Establish an ordering on alignments.
79 * Prefer aligned reads over unaligned reads
80 * Within the space of aligned reads
81 * prefer splice-free reads over splice reads, and
82 * indel-free reads over indel reads.
83 * If a read can either be indel-free or splice-free,
84 * prefer the indel-free alignment
85 */
86 bool AlignStatus::operator<(const AlignStatus& rhs) const
87 {
88 if (rhs._aligned != _aligned) return rhs._aligned;
89 if (rhs._edit_dist!=_edit_dist)
90 return rhs._edit_dist < _edit_dist;
91 //int lhs_value = _aligned ? 1 : 0;
92 int lhs_value = _indelFreeAlignment ? 4 : 0;
93 lhs_value += _unannotatedSpliceFreeAlignment ? 2 : 0;
94
95 //int rhs_value = rhs._aligned ? 1 : 0;
96 int rhs_value = rhs._indelFreeAlignment ? 4 : 0;
97 rhs_value += rhs._unannotatedSpliceFreeAlignment ? 2 : 0;
98 return lhs_value < rhs_value;
99 }
100
101 /**
102 * Alignments are only equal if all fields are identical.
103 */
104 bool AlignStatus::operator==(const AlignStatus& rhs) const
105 {
106 return ((_aligned == rhs._aligned) && (rhs._edit_dist ==_edit_dist) &&
107 (_indelFreeAlignment == rhs._indelFreeAlignment) &&
108 (_unannotatedSpliceFreeAlignment == rhs._unannotatedSpliceFreeAlignment));
109 }
110
111 bool AlignStatus::operator!=(const AlignStatus& rhs) const
112 {
113 return !((*this) == rhs);
114 }