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 (8 years, 4 months ago) by gpertea
File size: 3068 byte(s)
Log Message:
massive update with Daehwan's work

Line User Rev File contents
1 gpertea 29 /*
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 gpertea 135 _edit_dist = 0xFF;
33 gpertea 154 _fusionFreeAlignment = 0;
34 gpertea 29 }
35    
36     /**
37     * Parse the cigar string of a BowtieHit in order to determine the alignment status.
38     */
39 gpertea 154 AlignStatus::AlignStatus(const BowtieHit& bh, const JunctionSet& gtf_junctions){
40 gpertea 29 const vector<CigarOp>& cigar = bh.cigar();
41     _aligned = cigar.size() > 0;
42     _indelFreeAlignment = true;
43 gpertea 154 _fusionFreeAlignment = true;
44 gpertea 29 _unannotatedSpliceFreeAlignment = true;
45 gpertea 135 _edit_dist = bh.edit_dist();
46 gpertea 29 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 gpertea 154 case FUSION_FF:
73     case FUSION_FR:
74     case FUSION_RF:
75     case FUSION_RR:
76     _fusionFreeAlignment = false;
77     break;
78 gpertea 29 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 gpertea 135 if (rhs._aligned != _aligned) return rhs._aligned;
96     if (rhs._edit_dist!=_edit_dist)
97     return rhs._edit_dist < _edit_dist;
98 gpertea 154
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 gpertea 29
104 gpertea 154 // 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 gpertea 29 }
111    
112     /**
113     * Alignments are only equal if all fields are identical.
114     */
115     bool AlignStatus::operator==(const AlignStatus& rhs) const
116     {
117 gpertea 135 return ((_aligned == rhs._aligned) && (rhs._edit_dist ==_edit_dist) &&
118 gpertea 29 (_indelFreeAlignment == rhs._indelFreeAlignment) &&
119 gpertea 154 (_unannotatedSpliceFreeAlignment == rhs._unannotatedSpliceFreeAlignment) &&
120     (_fusionFreeAlignment == rhs._fusionFreeAlignment));
121 gpertea 29 }
122    
123     bool AlignStatus::operator!=(const AlignStatus& rhs) const
124     {
125     return !((*this) == rhs);
126     }