ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/tophat_cpp/deletions.cpp
Revision: 154
Committed: Tue Jan 24 02:29:21 2012 UTC (7 years, 8 months ago) by gpertea
File size: 4173 byte(s)
Log Message:
massive update with Daehwan's work

Line File contents
1 /*
2 * deletions.cpp
3 * TopHat
4 *
5 * Created by Ryan Kelley on 10/09/2010.
6 *
7 */
8
9 #ifdef HAVE_CONFIG_H
10 #include <config.h>
11 #endif
12
13 #include <cassert>
14 #include "common.h"
15 #include "deletions.h"
16
17
18
19
20
21 /*
22 * Print deletions in BED format
23 * As per the BED-standard (http://genome.ucsc.edu/FAQ/FAQformat)
24 * -The coordinates should be 0-based
25 * -The chromEnd field should not contain the actual feature
26 * -The name will be "-"
27 * -The score will be count of supporting reads (max of 1,000)
28 *
29 * chromStart refers to the position of the first deleted based
30 * <chrom>\t<left>\t<right>\t-\t<read count>\n
31 * @param deletions_out The output file
32 * @param deletions Maps from deletions to number of supporting reads
33 * @pram ref_sequences The table of reference sequences
34 *
35 */
36 void print_deletions(FILE* deletions_out, const DeletionSet& deletions, RefSequenceTable& ref_sequences){
37 fprintf(deletions_out, "track name=deletions description=\"TopHat deletions\"\n");
38 for(DeletionSet::const_iterator i = deletions.begin(); i != deletions.end(); ++i){
39 fprintf(deletions_out, "%s\t%d\t%d\t%s\t%d\n",
40 ref_sequences.get_name(i->first.refid),
41 i->first.left + 1,
42 i->first.right,
43 "-",
44 i->second);
45 }
46 }
47
48 /**
49 * Add deletions from an alignment to an DeletionSet.
50 * This will look for deletion in the alignment specified by bh. If the
51 * deletion is already in deletions, it will updated the count. Otherwise,
52 * it will add the deletion to the set and initialize the count to 1.
53 * @param bh The bowtie hit to be used to specify alignment infromation.
54 * @param deletions The DeletionSet that will be updated with the deletion information from teh alignment.
55 */
56 void deletions_from_alignment(const BowtieHit& bh, DeletionSet& deletions){
57 vector<Deletion> new_deletions;
58 deletions_from_spliced_hit(bh, new_deletions);
59
60 for(size_t i = 0; i < new_deletions.size(); ++i){
61 Deletion deletion = new_deletions[i];
62 DeletionSet::iterator itr = deletions.find(deletion);
63 if (itr != deletions.end()){
64 itr->second += 1;
65 }
66 else{
67 deletions[deletion] = 1;
68 }
69 }
70 return;
71 }
72
73
74
75
76 /**
77 * Extract a list of deletions from a bowtie hit.
78 * Given a bowtie hit, extract a vector of deletions.
79 * @param bh The bowtie hit to use for alignment information.
80 * @param insertions Used to store the resultant vector of deletions.
81 */
82 void deletions_from_spliced_hit(const BowtieHit& bh, vector<Deletion>& deletions){
83 const vector<CigarOp>& cigar = bh.cigar();
84 unsigned int positionInGenome = bh.left();
85 unsigned int positionInRead = 0;
86
87 bool bSawFusion = false;
88 for(size_t c = 0; c < cigar.size(); ++c){
89 Junction deletion;
90 switch(cigar[c].opcode){
91 case REF_SKIP:
92 positionInGenome += cigar[c].length;
93 break;
94 case rEF_SKIP:
95 positionInGenome -= cigar[c].length;
96 break;
97 case MATCH:
98 case mATCH:
99 if (cigar[c].opcode == MATCH)
100 positionInGenome += cigar[c].length;
101 else
102 positionInGenome -= cigar[c].length;
103 positionInRead += cigar[c].length;
104 break;
105 case DEL:
106 case dEL:
107 if (bSawFusion)
108 deletion.refid = bh.ref_id2();
109 else
110 deletion.refid = bh.ref_id();
111
112 if (cigar[c].opcode == DEL)
113 {
114 deletion.left = positionInGenome - 1;
115 deletion.right = positionInGenome + cigar[c].length;
116 }
117 else
118 {
119 deletion.left = positionInGenome - cigar[c].length;
120 deletion.right = positionInGenome + 1;
121 }
122
123 deletions.push_back(deletion);
124 positionInGenome += cigar[c].length;
125 break;
126 case INS:
127 case iNS:
128 positionInRead += cigar[c].length;
129 break;
130 case FUSION_FF:
131 case FUSION_FR:
132 case FUSION_RF:
133 bSawFusion = true;
134 positionInGenome = cigar[c].length;
135 break;
136 default:
137 break;
138 }
139 }
140 return;
141 }
142
143 void merge_with(DeletionSet& deletions, const DeletionSet& other)
144 {
145 for (DeletionSet::const_iterator deletion = other.begin(); deletion != other.end(); ++deletion)
146 {
147 DeletionSet::iterator itr = deletions.find(deletion->first);
148 if (itr != deletions.end())
149 {
150 itr->second += deletion->second;
151 }
152 else
153 {
154 deletions[deletion->first] = deletion->second;
155 }
156 }
157 }