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 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 29 }
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 gpertea 135 _edit_dist = bh.edit_dist();
45 gpertea 29 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 gpertea 135 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 gpertea 29 lhs_value += _unannotatedSpliceFreeAlignment ? 2 : 0;
94    
95 gpertea 135 //int rhs_value = rhs._aligned ? 1 : 0;
96     int rhs_value = rhs._indelFreeAlignment ? 4 : 0;
97 gpertea 29 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 gpertea 135 return ((_aligned == rhs._aligned) && (rhs._edit_dist ==_edit_dist) &&
107 gpertea 29 (_indelFreeAlignment == rhs._indelFreeAlignment) &&
108     (_unannotatedSpliceFreeAlignment == rhs._unannotatedSpliceFreeAlignment));
109     }
110    
111     bool AlignStatus::operator!=(const AlignStatus& rhs) const
112     {
113     return !((*this) == rhs);
114     }