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());
94 uint16_t excludeFlag = 0x0704,
95 uint16_t includeFlag = 0);
119 FUNC_CLASS myAnalyzeFuncPtr;
122 void addAlignmentPosition(
int refPosition,
SamRecord& record);
125 virtual void flushPileup(
int refID,
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)
279 processAlignment(record);
295 template <
class PILEUP_TYPE,
class FUNC_CLASS>
303 flushPileup(refID, refPosition);
310 addAlignmentPosition(refPosition, record);
315 template <
class PILEUP_TYPE,
class FUNC_CLASS>
326 flushPileup(refID, refPosition);
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))
375 flushPileup(pileupHead+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;
433 flushPileup(position);
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]);
Create/Access/Modify/Load Genome Sequences stored as binary mapped files.
Class to perform a pileup of all reads by position, assuming the reads are coordinate sorted.
virtual void processAlignment(SamRecord &record)
Add an alignment to the pileup.
virtual ~Pileup()
Destructor.
virtual int processFile(const std::string &fileName, uint16_t excludeFlag=0x0704, uint16_t includeFlag=0)
Performs a pileup on the specified file.
Pileup(int window, const std::string &refSeqFileName, const FUNC_CLASS &fp=FUNC_CLASS())
Perform pileup with a reference and a specified window size.
void flushPileup()
Done processing, flush every position that is currently being stored in the pileup.
Pileup(const FUNC_CLASS &fp=FUNC_CLASS())
Constructor using the default maximum number of bases a read spans.
Pileup(const std::string &refSeqFileName, const FUNC_CLASS &fp=FUNC_CLASS())
Perform pileup with a reference.
Pileup(int window, const FUNC_CLASS &fp=FUNC_CLASS())
Constructor that sets the maximum number of bases a read spans.
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...
Store refID/position, but does not store values < 0.
bool hasPosition(int refID, int refPosition)
Return whether or not this list contains the specified reference ID and position (negative values wil...
Allows the user to easily read/write a SAM/BAM file.
bool ReadHeader(SamFileHeader &header)
Reads the header section from the file and stores it in the passed in header.
const char * GetStatusMessage()
Get the Status Message of the last call that sets status.
@ COORDINATE
file is sorted by coordinate.
bool ReadRecord(SamFileHeader &header, SamRecord &record)
Reads the next record from the file & stores it in the passed in record.
void SetReference(GenomeSequence *reference)
Sets the reference to the specified genome sequence object.
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...
SamStatus::Status GetStatus()
Get the Status of the last call that sets status.
void setSortedValidation(SortedType sortType)
Set the flag to validate that the file is sorted as it is read/written.
Class providing an easy to use interface to get/set/operate on the fields in a SAM/BAM record.
int32_t get1BasedPosition()
Get the 1-based(SAM) leftmost position (POS) of the record.
int32_t getReferenceID()
Get the reference sequence id of the record (BAM format rid).
uint16_t getFlag()
Get the flag (FLAG).
int32_t get0BasedAlignmentEnd()
Returns the 0-based inclusive rightmost position of the clipped sequence.
int32_t get0BasedPosition()
Get the 0-based(BAM) leftmost position of the record.
const char * getCigar()
Returns the SAM formatted CIGAR string.
const char * getReadName()
Returns the SAM formatted Read Name (QNAME).
@ NO_MORE_RECS
NO_MORE_RECS: failed to read a record since there are no more to read either in the file or section i...