ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/tophat_cpp/insertions.cpp
(Generate patch)
# Line 110 | Line 110
110          unsigned int positionInGenome = bh.left();
111          unsigned int positionInRead = 0;
112  
113 +        bool bSawFusion = false;
114          for(size_t c = 0; c < cigar.size(); ++c){
115 <                Insertion insertion;
116 <                switch(cigar[c].opcode){
117 <                        case REF_SKIP:
118 <                                positionInGenome += cigar[c].length;
119 <                                break;
120 <                        case MATCH:
121 <                                positionInGenome += cigar[c].length;
122 <                                positionInRead += cigar[c].length;
123 <                                break;
124 <                        case DEL:
125 <                                positionInGenome += cigar[c].length;
126 <                                break;
127 <                        case INS:
128 <                                /*
129 <                                 * Note that the reported position in the genome from the SAM
130 <                                 * alignment is 1-based, since the insertion object is expecting
131 <                                 * a 0-based co-ordinate, we need to subtract 1
132 <                                 */
133 <                                insertion.refid = bh.ref_id();
134 <                                insertion.left = positionInGenome-1;
135 <                                insertion.sequence = bh.seq().substr(positionInRead, cigar[c].length);
115 >          Insertion insertion;
116 >          switch(cigar[c].opcode){
117 >          case REF_SKIP:
118 >            positionInGenome += cigar[c].length;
119 >            break;
120 >          case rEF_SKIP:
121 >            positionInGenome -= cigar[c].length;
122 >            break;
123 >          case MATCH:
124 >          case mATCH:
125 >            if (cigar[c].opcode == MATCH)
126 >              positionInGenome += cigar[c].length;
127 >            else
128 >              positionInGenome -= cigar[c].length;
129 >            positionInRead += cigar[c].length;
130 >            break;
131 >          case DEL:
132 >            positionInGenome += cigar[c].length;
133 >            break;
134 >          case dEL:
135 >            positionInGenome -= cigar[c].length;
136 >            break;
137 >          case INS:
138 >          case iNS:
139 >            /*
140 >             * Note that the reported position in the genome from the SAM
141 >             * alignment is 1-based, since the insertion object is expecting
142 >             * a 0-based co-ordinate, we need to subtract 1
143 >             */
144 >            if (bSawFusion)
145 >              insertion.refid = bh.ref_id2();
146 >            else
147 >              insertion.refid = bh.ref_id();
148  
149 <                                insertions.push_back(insertion);
150 <                                positionInRead += cigar[c].length;
151 <                                break;
152 <                        default:
153 <                                break;
149 >            if (cigar[c].opcode == INS)
150 >              insertion.left = positionInGenome - 1;
151 >            else
152 >              insertion.left = positionInGenome + 1;
153 >              
154 >            insertion.sequence = bh.seq().substr(positionInRead, cigar[c].length);
155 >            
156 >            insertions.push_back(insertion);
157 >            positionInRead += cigar[c].length;
158 >            break;
159 >          case FUSION_FF:
160 >          case FUSION_FR:
161 >          case FUSION_RF:
162 >            bSawFusion = true;
163 >            positionInGenome = cigar[c].length;
164 >            break;
165 >          default:
166 >            break;
167                  }      
168          }      
169          return;
170   }
171  
172 <
173 <
172 > void merge_with(InsertionSet& insertions, const InsertionSet& other)
173 > {
174 >  for (InsertionSet::const_iterator insertion = other.begin(); insertion != other.end(); ++insertion)
175 >    {
176 >      InsertionSet::iterator itr = insertions.find(insertion->first);
177 >      if (itr != insertions.end())
178 >        {
179 >          itr->second += insertion->second;
180 >        }
181 >      else
182 >        {
183 >          insertions[insertion->first] = insertion->second;
184 >        }
185 >    }
186 > }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines