00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #ifndef __PILEUP_H__
00019 #define __PILEUP_H__
00020
00021 template <class PILEUP_TYPE>
00022 class defaultPileup
00023 {
00024 public:
00025 bool operator() (PILEUP_TYPE& element)
00026 {
00027 element.analyze();
00028 return(true);
00029 }
00030 void analyze(PILEUP_TYPE element)
00031 {
00032 element.analyze();
00033 }
00034 };
00035
00036 template <class PILEUP_TYPE>
00037 void defaultPileupAnalyze(PILEUP_TYPE& element)
00038 {
00039 element.analyze();
00040 }
00041
00042
00043 template <class PILEUP_TYPE,
00044 class FUNC_CLASS = defaultPileup<PILEUP_TYPE> >
00045 class Pileup
00046 {
00047 public:
00048 Pileup(const FUNC_CLASS& fp = FUNC_CLASS());
00049
00050 Pileup(int window, const FUNC_CLASS& fp = FUNC_CLASS());
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063 virtual int processFile(const std::string& fileName,
00064 uint16_t excludeFlag = 0x0704,
00065 uint16_t includeFlag = 0);
00066
00067
00068
00069
00070 virtual void processAlignment(SamRecord& record);
00071
00072
00073 void flushPileup();
00074
00075 protected:
00076 FUNC_CLASS myAnalyzeFuncPtr;
00077
00078 private:
00079
00080
00081 void addAlignmentPosition(int refPosition, SamRecord& record);
00082
00083
00084 void flushPileup(int refID, int refPosition);
00085 void flushPileup(int refPosition);
00086
00087
00088
00089
00090 int pileupPosition(int refPosition);
00091
00092 virtual void resetElement(PILEUP_TYPE& element, int position);
00093 virtual void addElement(PILEUP_TYPE& element, SamRecord& record);
00094 virtual void analyzeElement(PILEUP_TYPE& element);
00095
00096
00097 std::vector<PILEUP_TYPE> myElements;
00098
00099 int pileupStart;
00100 int pileupHead;
00101 int pileupTail;
00102 int pileupWindow;
00103
00104 int myCurrentRefID;
00105 };
00106
00107
00108 template <class PILEUP_TYPE, class FUNC_CLASS>
00109 Pileup<PILEUP_TYPE, FUNC_CLASS>::Pileup(const FUNC_CLASS& fp)
00110 : myAnalyzeFuncPtr(fp),
00111 myElements(),
00112 pileupStart(0),
00113 pileupHead(0),
00114 pileupTail(0),
00115 pileupWindow(1024),
00116 myCurrentRefID(-2)
00117 {
00118
00119 myElements.resize(pileupWindow, PILEUP_TYPE());
00120 }
00121
00122
00123 template <class PILEUP_TYPE, class FUNC_CLASS>
00124 Pileup<PILEUP_TYPE, FUNC_CLASS>::Pileup(int window, const FUNC_CLASS& fp)
00125 : myAnalyzeFuncPtr(fp),
00126 myElements(),
00127 pileupStart(0),
00128 pileupHead(0),
00129 pileupTail(0),
00130 pileupWindow(window),
00131 myCurrentRefID(-2)
00132 {
00133
00134 myElements.resize(window, PILEUP_TYPE());
00135 }
00136
00137
00138 template <class PILEUP_TYPE, class FUNC_CLASS>
00139 int Pileup<PILEUP_TYPE, FUNC_CLASS>::processFile(const std::string& fileName,
00140 uint16_t excludeFlag,
00141 uint16_t includeFlag)
00142 {
00143 SamFile samIn;
00144 SamFileHeader header;
00145 SamRecord record;
00146
00147 if(!samIn.OpenForRead(fileName.c_str()))
00148 {
00149 fprintf(stderr, "%s\n", samIn.GetStatusMessage());
00150 return(samIn.GetStatus());
00151 }
00152
00153 if(!samIn.ReadHeader(header))
00154 {
00155 fprintf(stderr, "%s\n", samIn.GetStatusMessage());
00156 return(samIn.GetStatus());
00157 }
00158
00159
00160 samIn.setSortedValidation(SamFile::COORDINATE);
00161
00162
00163 while (samIn.ReadRecord(header, record))
00164 {
00165 uint16_t flag = record.getFlag();
00166 if(flag & excludeFlag)
00167 {
00168
00169
00170 continue;
00171 }
00172 if((flag & includeFlag) != includeFlag)
00173 {
00174
00175
00176 continue;
00177 }
00178 processAlignment(record);
00179 }
00180
00181 flushPileup();
00182
00183 int returnValue = 0;
00184 if(samIn.GetStatus() != SamStatus::NO_MORE_RECS)
00185 {
00186
00187 fprintf(stderr, "%s\n", samIn.GetStatusMessage());
00188 returnValue = samIn.GetStatus();
00189 }
00190 return(returnValue);
00191 }
00192
00193
00194 template <class PILEUP_TYPE, class FUNC_CLASS>
00195 void Pileup<PILEUP_TYPE, FUNC_CLASS>::processAlignment(SamRecord& record)
00196 {
00197 int refPosition = record.get0BasedPosition();
00198 int refID = record.getReferenceID();
00199
00200
00201
00202 flushPileup(refID, refPosition);
00203
00204
00205
00206
00207 for(; refPosition <= record.get0BasedAlignmentEnd(); ++refPosition)
00208 {
00209 addAlignmentPosition(refPosition, record);
00210 }
00211 }
00212
00213
00214 template <class PILEUP_TYPE, class FUNC_CLASS>
00215 void Pileup<PILEUP_TYPE, FUNC_CLASS>::flushPileup()
00216 {
00217 while (pileupHead <= pileupTail)
00218 {
00219 flushPileup(pileupHead+1);
00220 }
00221 pileupStart = pileupHead = pileupTail = 0;
00222 }
00223
00224
00225
00226 template <class PILEUP_TYPE, class FUNC_CLASS>
00227 void Pileup<PILEUP_TYPE, FUNC_CLASS>::addAlignmentPosition(int refPosition,
00228 SamRecord& record)
00229 {
00230 int offset = pileupPosition(refPosition);
00231
00232 if((offset < 0) || (offset >= pileupWindow))
00233 {
00234 std::cerr << "Pileup Buffer Overflow: position = " << refPosition
00235 << "; refID = " << record.getReferenceID()
00236 << "; recStartPos = " << record.get1BasedPosition()
00237 << "; pileupStart = " << pileupStart
00238 << "; pileupHead = " << pileupHead
00239 << "; pileupTail = " << pileupTail;
00240 }
00241
00242 addElement(myElements[offset], record);
00243 }
00244
00245
00246 template <class PILEUP_TYPE, class FUNC_CLASS>
00247 void Pileup<PILEUP_TYPE, FUNC_CLASS>::flushPileup(int refID, int position)
00248 {
00249
00250 if(refID != myCurrentRefID)
00251 {
00252
00253 flushPileup();
00254 myCurrentRefID = refID;
00255
00256
00257 pileupStart = pileupHead = position;
00258 }
00259 else
00260 {
00261
00262 flushPileup(position);
00263 }
00264 }
00265
00266
00267 template <class PILEUP_TYPE, class FUNC_CLASS>
00268 void Pileup<PILEUP_TYPE, FUNC_CLASS>::flushPileup(int position)
00269 {
00270
00271 while((pileupHead < position) && (pileupHead <= pileupTail))
00272 {
00273 analyzeElement(myElements[pileupHead - pileupStart]);
00274
00275 pileupHead++;
00276
00277 if(pileupHead - pileupStart >= pileupWindow)
00278 pileupStart += pileupWindow;
00279 }
00280
00281
00282
00283 if(pileupHead != position)
00284 {
00285 pileupHead = pileupStart = position;
00286 }
00287 }
00288
00289
00290
00291
00292
00293 template <class PILEUP_TYPE, class FUNC_CLASS>
00294 int Pileup<PILEUP_TYPE, FUNC_CLASS>::pileupPosition(int position)
00295 {
00296 if((position < pileupHead) || (position > (pileupHead + pileupWindow)))
00297 {
00298 return -1;
00299 error("Pileup Buffer Overflow");
00300 }
00301
00302 int offset = position - pileupStart;
00303
00304 if(offset >= pileupWindow)
00305 {
00306 offset -= pileupWindow;
00307 }
00308
00309 if(position > pileupTail)
00310 {
00311
00312
00313 resetElement(myElements[offset], position);
00314 pileupTail = position;
00315 }
00316
00317 return(offset);
00318 }
00319
00320
00321 template <class PILEUP_TYPE, class FUNC_CLASS>
00322 void Pileup<PILEUP_TYPE, FUNC_CLASS>::resetElement(PILEUP_TYPE& element,
00323 int position)
00324 {
00325 element.reset(position);
00326 }
00327
00328
00329 template <class PILEUP_TYPE, class FUNC_CLASS>
00330 void Pileup<PILEUP_TYPE, FUNC_CLASS>::addElement(PILEUP_TYPE& element,
00331 SamRecord& record)
00332 {
00333 element.addEntry(record);
00334 }
00335
00336
00337 template <class PILEUP_TYPE, class FUNC_CLASS>
00338 void Pileup<PILEUP_TYPE, FUNC_CLASS>::analyzeElement(PILEUP_TYPE& element)
00339 {
00340 myAnalyzeFuncPtr(element);
00341 }
00342
00343
00344
00345
00346 #endif