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 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     }
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     }