29 static const int DEFAULT_WINDOW_SIZE = 1024;
32 template <
class PILEUP_TYPE>
36 bool operator() (PILEUP_TYPE& element)
41 void analyze(PILEUP_TYPE element)
47 template <
class PILEUP_TYPE>
48 void defaultPileupAnalyze(PILEUP_TYPE& element)
56 template <
class PILEUP_TYPE,
62 Pileup(
const FUNC_CLASS& fp = FUNC_CLASS());
68 Pileup(
int window,
const FUNC_CLASS& fp = FUNC_CLASS());
71 Pileup(
const std::string& refSeqFileName,
72 const FUNC_CLASS& fp = FUNC_CLASS());
75 Pileup(
int window,
const std::string& refSeqFileName,
76 const FUNC_CLASS& fp = FUNC_CLASS());
93 virtual int processFile(
const std::string& fileName,
94 uint16_t excludeFlag = 0x0704,
95 uint16_t includeFlag = 0);
98 virtual void processAlignment(
SamRecord& record);
109 virtual void processAlignmentRegion(
SamRecord& record,
119 FUNC_CLASS myAnalyzeFuncPtr;
122 void addAlignmentPosition(
int refPosition,
SamRecord& record);
125 virtual void flushPileup(
int refID,
int refPosition);
126 void flushPileup(
int refPosition);
131 int pileupPosition(
int refPosition);
133 virtual void resetElement(PILEUP_TYPE& element,
int position);
134 virtual void addElement(PILEUP_TYPE& element,
SamRecord& record);
135 virtual void analyzeElement(PILEUP_TYPE& element);
136 virtual void analyzeHead();
138 std::vector<PILEUP_TYPE> myElements;
151 template <
class PILEUP_TYPE,
class FUNC_CLASS>
153 : myAnalyzeFuncPtr(fp),
163 myElements.resize(pileupWindow);
167 template <
class PILEUP_TYPE,
class FUNC_CLASS>
169 : myAnalyzeFuncPtr(fp),
174 pileupWindow(window),
179 myElements.resize(window);
183 template <
class PILEUP_TYPE,
class FUNC_CLASS>
185 : myAnalyzeFuncPtr(fp),
197 myElements.resize(pileupWindow);
199 PILEUP_TYPE::setReference(myRefPtr);
203 template <
class PILEUP_TYPE,
class FUNC_CLASS>
205 : myAnalyzeFuncPtr(fp),
210 pileupWindow(window),
217 myElements.resize(window);
219 PILEUP_TYPE::setReference(myRefPtr);
223 template <
class PILEUP_TYPE,
class FUNC_CLASS>
234 template <
class PILEUP_TYPE,
class FUNC_CLASS>
236 uint16_t excludeFlag,
237 uint16_t includeFlag)
266 uint16_t flag = record.
getFlag();
267 if(flag & excludeFlag)
273 if((flag & includeFlag) != includeFlag)
295 template <
class PILEUP_TYPE,
class FUNC_CLASS>
310 addAlignmentPosition(refPosition, record);
315 template <
class PILEUP_TYPE,
class FUNC_CLASS>
330 if(startPos > refPosition)
332 refPosition = startPos;
343 if((endPos != -1) && (refPosition >= endPos))
350 if(excludeList != NULL)
361 addAlignmentPosition(refPosition, record);
367 template <
class PILEUP_TYPE,
class FUNC_CLASS>
373 while ((pileupHead <= pileupTail) && (pileupTail != -1))
377 pileupStart = pileupHead = 0;
383 template <
class PILEUP_TYPE,
class FUNC_CLASS>
389 offset = pileupPosition(refPosition);
391 catch(std::runtime_error& err)
393 const char* overflowErr =
"Overflow on the pileup buffer:";
394 String errorMessage = err.what();
395 if(strncmp(err.what(), overflowErr, strlen(overflowErr)) == 0)
398 errorMessage +=
"\n\tPileup Buffer Overflow: recordName = ";
400 errorMessage +=
"; Cigar = ";
403 throw std::runtime_error(errorMessage.c_str());
406 if((offset < 0) || (offset >= pileupWindow))
408 std::cerr <<
"Pileup Buffer Overflow: position = " << refPosition
411 <<
"; pileupStart = " << pileupStart
412 <<
"; pileupHead = " << pileupHead
413 <<
"; pileupTail = " << pileupTail;
416 addElement(myElements[offset], record);
420 template <
class PILEUP_TYPE,
class FUNC_CLASS>
424 if(refID != myCurrentRefID)
428 myCurrentRefID = refID;
438 template <
class PILEUP_TYPE,
class FUNC_CLASS>
443 while((pileupHead < position) && (pileupHead <= pileupTail))
449 if(pileupHead - pileupStart >= pileupWindow)
450 pileupStart += pileupWindow;
453 if(pileupHead > pileupTail)
456 pileupHead = pileupStart = 0;
465 template <
class PILEUP_TYPE,
class FUNC_CLASS>
471 pileupStart = pileupHead = position;
474 resetElement(myElements[0], position);
475 pileupTail = position;
480 if((position < pileupHead) || (position > (pileupHead + pileupWindow)))
483 "Overflow on the pileup buffer: specifiedPosition = ";
484 errorMessage += position;
485 errorMessage +=
", pileup buffer start position: ";
486 errorMessage += pileupHead;
487 errorMessage +=
", pileup buffer end position: ";
488 errorMessage += pileupHead + pileupWindow;
490 throw std::runtime_error(errorMessage.c_str());
494 int offset = position - pileupStart;
496 if(offset >= pileupWindow)
498 offset -= pileupWindow;
503 while(position > pileupTail)
510 offset = pileupTail - pileupStart;
511 if(offset >= pileupWindow)
513 offset -= pileupWindow;
518 resetElement(myElements[offset], pileupTail);
525 template <
class PILEUP_TYPE,
class FUNC_CLASS>
529 element.reset(position);
533 template <
class PILEUP_TYPE,
class FUNC_CLASS>
537 element.addEntry(record);
541 template <
class PILEUP_TYPE,
class FUNC_CLASS>
544 myAnalyzeFuncPtr(element);
548 template <
class PILEUP_TYPE,
class FUNC_CLASS>
551 myAnalyzeFuncPtr(myElements[pileupHead - pileupStart]);
int32_t get0BasedAlignmentEnd()
Returns the 0-based inclusive rightmost position of the clipped sequence.
bool ReadHeader(SamFileHeader &header)
Reads the header section from the file and stores it in the passed in header.
void SetReference(GenomeSequence *reference)
Sets the reference to the specified genome sequence object.
virtual ~Pileup()
Destructor.
const char * getReadName()
Returns the SAM formatted Read Name (QNAME).
int32_t getReferenceID()
Get the reference sequence id of the record (BAM format rid).
void setSortedValidation(SortedType sortType)
Set the flag to validate that the file is sorted as it is read/written.
NO_MORE_RECS: failed to read a record since there are no more to read either in the file or section i...
file is sorted by coordinate.
void flushPileup()
Done processing, flush every position that is currently being stored in the pileup.
int32_t get1BasedPosition()
Get the 1-based(SAM) leftmost position (POS) of the record.
Allows the user to easily read/write a SAM/BAM file.
SamStatus::Status GetStatus()
Get the Status of the last call that sets status.
virtual void processAlignmentRegion(SamRecord &record, int startPos, int endPos, PosList *excludeList=NULL)
Add only positions that fall within the specified region of the alignment to the pileup and outside o...
virtual void processAlignment(SamRecord &record)
Add an alignment to the pileup.
Pileup(const FUNC_CLASS &fp=FUNC_CLASS())
Constructor using the default maximum number of bases a read spans.
const char * GetStatusMessage()
Get the Status Message of the last call that sets status.
bool ReadRecord(SamFileHeader &header, SamRecord &record)
Reads the next record from the file & stores it in the passed in record.
Create/Access/Modify/Load Genome Sequences stored as binary mapped files.
Store refID/position, but does not store values < 0.
Class providing an easy to use interface to get/set/operate on the fields in a SAM/BAM record...
bool OpenForRead(const char *filename, SamFileHeader *header=NULL)
Open a sam/bam file for reading with the specified filename, determing the type of file and SAM/BAM b...
uint16_t getFlag()
Get the flag (FLAG).
int32_t get0BasedPosition()
Get the 0-based(BAM) leftmost position of the record.
virtual int processFile(const std::string &fileName, uint16_t excludeFlag=0x0704, uint16_t includeFlag=0)
Performs a pileup on the specified file.
bool hasPosition(int refID, int refPosition)
Return whether or not this list contains the specified reference ID and position (negative values wil...
const char * getCigar()
Returns the SAM formatted CIGAR string.
Class to perform a pileup of all reads by position, assuming the reads are coordinate sorted...