CigarHelper.cpp
00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include "CigarHelper.h"
00021
00022
00023 int32_t CigarHelper::softClipBeginByRefPos(SamRecord& record,
00024 int32_t refPosition0Based,
00025 CigarRoller& newCigar,
00026 int32_t &new0BasedPosition)
00027 {
00028 newCigar.clear();
00029 Cigar* cigar = record.getCigarInfo();
00030 if(cigar == NULL)
00031 {
00032
00033 ErrorHandler::handleError("Soft clipping, but failed to read the cigar");
00034 return(NO_CLIP);
00035 }
00036
00037
00038 if((cigar->size() == 0) || (record.get0BasedPosition() == -1))
00039 {
00040 return(NO_CLIP);
00041 }
00042
00043
00044
00045 if(refPosition0Based < record.get0BasedPosition())
00046 {
00047
00048 newCigar.Set(record.getCigar());
00049 return(NO_CLIP);
00050 }
00051
00052
00053
00054 int32_t readClipPosition = 0;
00055 bool clipWritten = false;
00056 new0BasedPosition = record.get0BasedPosition();
00057 for(int i = 0; i < cigar->size(); i++)
00058 {
00059 const Cigar::CigarOperator* op = &(cigar->getOperator(i));
00060
00061 if(clipWritten)
00062 {
00063
00064 newCigar += *op;
00065
00066 continue;
00067 }
00068
00069
00070
00071
00072
00073
00074 if(Cigar::foundInReference(*op))
00075 {
00076
00077
00078
00079
00080 new0BasedPosition += op->count;
00081
00082
00083
00084 if(Cigar::foundInQuery(*op))
00085 {
00086
00087
00088
00089 uint32_t numKeep = 0;
00090
00091 if(refPosition0Based < new0BasedPosition)
00092 {
00093
00094 numKeep = new0BasedPosition - refPosition0Based - 1;
00095
00096 if(numKeep > op->count)
00097 {
00098
00099
00100
00101 numKeep = op->count;
00102 }
00103 }
00104
00105
00106
00107 readClipPosition += (op->count - numKeep);
00108
00109
00110
00111
00112 if(numKeep > 0)
00113 {
00114 new0BasedPosition -= numKeep;
00115
00116 newCigar.Add(Cigar::softClip, readClipPosition);
00117
00118
00119
00120 newCigar.Add(op->operation, numKeep);
00121
00122
00123
00124 clipWritten = true;
00125 continue;
00126 }
00127 }
00128 }
00129 else
00130 {
00131
00132
00133 if(op->operation == Cigar::hardClip)
00134 {
00135
00136 if(i == 0)
00137 {
00138 newCigar += *op;
00139 }
00140
00141 else if(i == (cigar->size() - 1))
00142 {
00143
00144
00145 if(clipWritten == false)
00146 {
00147 newCigar.Add(Cigar::softClip, readClipPosition);
00148
00149
00150 new0BasedPosition = record.get0BasedPosition();
00151 clipWritten = true;
00152 }
00153
00154 newCigar += *op;
00155 }
00156 }
00157
00158 if(Cigar::foundInQuery(*op))
00159 {
00160
00161 readClipPosition += op->count;
00162 }
00163 }
00164 }
00165
00166
00167
00168
00169 if(clipWritten == false)
00170 {
00171 newCigar.Add(Cigar::softClip, readClipPosition);
00172
00173
00174 new0BasedPosition = record.get0BasedPosition();
00175 }
00176
00177
00178
00179 return(readClipPosition - 1);
00180 }
00181
00182
00183
00184 int32_t CigarHelper::softClipEndByRefPos(SamRecord& record,
00185 int32_t refPosition0Based,
00186 CigarRoller& newCigar)
00187 {
00188 newCigar.clear();
00189 Cigar* cigar = record.getCigarInfo();
00190 if(cigar == NULL)
00191 {
00192
00193 ErrorHandler::handleError("Soft clipping, but failed to read the cigar");
00194 return(NO_CLIP);
00195 }
00196
00197
00198 if((cigar->size() == 0) || (record.get0BasedPosition() == -1))
00199 {
00200 return(NO_CLIP);
00201 }
00202
00203
00204
00205 if(refPosition0Based > record.get0BasedAlignmentEnd())
00206 {
00207
00208 newCigar.Set(record.getCigar());
00209 return(NO_CLIP);
00210 }
00211
00212
00213
00214 int32_t currentRefPosition = record.get0BasedPosition();
00215 int32_t readClipPosition = 0;
00216 for(int i = 0; i < cigar->size(); i++)
00217 {
00218 const Cigar::CigarOperator* op = &(cigar->getOperator(i));
00219
00220
00221
00222 if(Cigar::foundInReference(*op))
00223 {
00224
00225
00226
00227 currentRefPosition += op->count;
00228 }
00229
00230
00231 if(refPosition0Based < currentRefPosition)
00232 {
00233
00234
00235 int32_t numKeep = 0;
00236 if(Cigar::foundInQuery(*op))
00237 {
00238 numKeep = op->count - (currentRefPosition - refPosition0Based);
00239 if(numKeep > 0)
00240 {
00241 newCigar.Add(op->operation, numKeep);
00242 readClipPosition += numKeep;
00243 }
00244 }
00245 else if(Cigar::isClip(*op))
00246 {
00247
00248 newCigar.Add(op->operation, op->count);
00249 }
00250 else
00251 {
00252
00253
00254
00255 }
00256
00257 break;
00258 }
00259 else if(refPosition0Based == currentRefPosition)
00260 {
00261 newCigar += *op;
00262 if(Cigar::foundInQuery(*op))
00263 {
00264 readClipPosition += op->count;
00265 }
00266 }
00267 else
00268 {
00269
00270
00271 newCigar += *op;
00272 if(Cigar::foundInQuery(*op))
00273 {
00274
00275 readClipPosition += op->count;
00276 }
00277 }
00278 }
00279
00280
00281
00282
00283
00284 for(int j = newCigar.size() - 1; j >= 0; j--)
00285 {
00286 const Cigar::CigarOperator* op = &(newCigar.getOperator(j));
00287 if(!Cigar::foundInQuery(*op) && !Cigar::isClip(*op))
00288 {
00289
00290 newCigar.Remove(j);
00291 }
00292 else if(Cigar::foundInQuery(*op) & Cigar::isClip(*op))
00293 {
00294
00295
00296
00297 readClipPosition -= op->count;
00298 newCigar.Remove(j);
00299 }
00300 else
00301 {
00302
00303 break;
00304 }
00305 }
00306
00307
00308 int32_t numSoftClips = record.getReadLength() - readClipPosition;
00309
00310
00311 newCigar.Add(Cigar::softClip, numSoftClips);
00312
00313
00314 if(cigar->size() != 0)
00315 {
00316 const Cigar::CigarOperator* lastOp =
00317 &(cigar->getOperator(cigar->size() - 1));
00318 if(lastOp->operation == Cigar::hardClip)
00319 {
00320 newCigar += *lastOp;
00321 }
00322 }
00323
00324 return(readClipPosition);
00325 }