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