Class to perform a pileup of all reads by position, assuming the reads are coordinate sorted. More...

#include <Pileup.h>

Collaboration diagram for Pileup< PILEUP_TYPE, FUNC_CLASS >:

Public Member Functions
	Pileup (const FUNC_CLASS &fp=FUNC_CLASS())
	Constructor using the default maximum number of bases a read spans.
	Pileup (int window, const FUNC_CLASS &fp=FUNC_CLASS())
	Constructor that sets the 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 std::string &refSeqFileName, const FUNC_CLASS &fp=FUNC_CLASS())
	Perform pileup with a reference and a specified window size.
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.
virtual void	processAlignment (SamRecord &record)
	Add an alignment to the pileup.
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 of the specified excluded positions.
void	flushPileup ()
	Done processing, flush every position that is currently being stored in the pileup.

Protected Member Functions
void	addAlignmentPosition (int refPosition, SamRecord &record)
virtual void	flushPileup (int refID, int refPosition)
void	flushPileup (int refPosition)
int	pileupPosition (int refPosition)
virtual void	resetElement (PILEUP_TYPE &element, int position)
virtual void	addElement (PILEUP_TYPE &element, SamRecord &record)
virtual void	analyzeElement (PILEUP_TYPE &element)
virtual void	analyzeHead ()

Protected Attributes
FUNC_CLASS	myAnalyzeFuncPtr
std::vector< PILEUP_TYPE >	myElements
int	pileupStart
int	pileupHead
int	pileupTail
int	pileupWindow
int	myCurrentRefID
GenomeSequence *	myRefPtr

Detailed Description

template<class PILEUP_TYPE, class FUNC_CLASS = defaultPileup<PILEUP_TYPE>>
class Pileup< PILEUP_TYPE, FUNC_CLASS >

Class to perform a pileup of all reads by position, assuming the reads are coordinate sorted.

Definition at line 58 of file Pileup.h.

Constructor & Destructor Documentation

◆ Pileup() [1/4]

template<class PILEUP_TYPE, class FUNC_CLASS>

Pileup< PILEUP_TYPE, FUNC_CLASS >::Pileup ( const FUNC_CLASS & fp = FUNC_CLASS() )

Constructor using the default maximum number of bases a read spans.

Definition at line 152 of file Pileup.h.

    : myAnalyzeFuncPtr(fp),
      myElements(),
      pileupStart(0),
      pileupHead(0),
      pileupTail(-1),
      pileupWindow(PileupHelper::DEFAULT_WINDOW_SIZE),
      myCurrentRefID(-2),
      myRefPtr(NULL)
{
    // Not using pointers since this is templated.
    myElements.resize(pileupWindow);
}

◆ Pileup() [2/4]

template<class PILEUP_TYPE, class FUNC_CLASS>

Pileup< PILEUP_TYPE, FUNC_CLASS >::Pileup	(	int	window,
		const FUNC_CLASS &	fp = FUNC_CLASS() )

Constructor that sets the maximum number of bases a read spans.

This is the "window" the length of the buffer that holds the pileups for each position until the read start has moved past the position.

Definition at line 168 of file Pileup.h.

    : myAnalyzeFuncPtr(fp),
      myElements(),
      pileupStart(0),
      pileupHead(0),
      pileupTail(-1),
      pileupWindow(window),
      myCurrentRefID(-2),
      myRefPtr(NULL)
{
    // Not using pointers since this is templated.
    myElements.resize(window);
}

◆ Pileup() [3/4]

template<class PILEUP_TYPE, class FUNC_CLASS>

Pileup< PILEUP_TYPE, FUNC_CLASS >::Pileup	(	const std::string &	refSeqFileName,
		const FUNC_CLASS &	fp = FUNC_CLASS() )

Perform pileup with a reference.

Definition at line 184 of file Pileup.h.

    : myAnalyzeFuncPtr(fp),
      myElements(),
      pileupStart(0),
      pileupHead(0),
      pileupTail(-1),
      pileupWindow(PileupHelper::DEFAULT_WINDOW_SIZE),
      myCurrentRefID(-2),
      myRefPtr(NULL)
{
    myRefPtr = new GenomeSequence(refSeqFileName.c_str());
 
    // Not using pointers since this is templated.
    myElements.resize(pileupWindow);
 
    PILEUP_TYPE::setReference(myRefPtr);
}

◆ Pileup() [4/4]

template<class PILEUP_TYPE, class FUNC_CLASS>

Pileup< PILEUP_TYPE, FUNC_CLASS >::Pileup	(	int	window,
		const std::string &	refSeqFileName,
		const FUNC_CLASS &	fp = FUNC_CLASS() )

Perform pileup with a reference and a specified window size.

Definition at line 204 of file Pileup.h.

    : myAnalyzeFuncPtr(fp),
      myElements(),
      pileupStart(0),
      pileupHead(0),
      pileupTail(-1),
      pileupWindow(window),
      myCurrentRefID(-2),
      myRefPtr(NULL)
{
    myRefPtr = new GenomeSequence(refSeqFileName.c_str());
 
    // Not using pointers since this is templated.
    myElements.resize(window);
 
    PILEUP_TYPE::setReference(myRefPtr);
}

◆ ~Pileup()

template<class PILEUP_TYPE, class FUNC_CLASS>

Pileup< PILEUP_TYPE, FUNC_CLASS >::~Pileup ( )

virtual

Destructor.

Definition at line 224 of file Pileup.h.

{
    flushPileup();
    if(myRefPtr != NULL)
    {
        delete myRefPtr;
        myRefPtr = NULL;
    }
}

References flushPileup().

Member Function Documentation

◆ addAlignmentPosition()

template<class PILEUP_TYPE, class FUNC_CLASS>

void Pileup< PILEUP_TYPE, FUNC_CLASS >::addAlignmentPosition	(	int	refPosition,
		SamRecord &	record )

protected

Definition at line 384 of file Pileup.h.

{
    int offset = 0;
    try{
        offset = pileupPosition(refPosition);
    }
    catch(std::runtime_error& err)
    {
        const char* overflowErr = "Overflow on the pileup buffer:";
        String errorMessage = err.what();
        if(strncmp(err.what(), overflowErr, strlen(overflowErr)) == 0)
        {
            
            errorMessage += "\n\tPileup Buffer Overflow: recordName = ";
            errorMessage += record.getReadName();
            errorMessage += "; Cigar = ";
            errorMessage += record.getCigar();
        }
        throw std::runtime_error(errorMessage.c_str());
    }
    
    if((offset < 0) || (offset >= pileupWindow))
    {
        std::cerr << "Pileup Buffer Overflow: position = " << refPosition
                  << "; refID = " << record.getReferenceID() 
                  << "; recStartPos = " << record.get1BasedPosition()
                  << "; pileupStart = " << pileupStart
                  << "; pileupHead = " << pileupHead
                  << "; pileupTail = " << pileupTail;
    }
 
    addElement(myElements[offset], record);
}

◆ addElement()

template<class PILEUP_TYPE, class FUNC_CLASS>

void Pileup< PILEUP_TYPE, FUNC_CLASS >::addElement	(	PILEUP_TYPE &	element,
		SamRecord &	record )

protectedvirtual

Definition at line 534 of file Pileup.h.

{
    element.addEntry(record);
}

◆ analyzeElement()

template<class PILEUP_TYPE, class FUNC_CLASS>

void Pileup< PILEUP_TYPE, FUNC_CLASS >::analyzeElement ( PILEUP_TYPE & element )

protectedvirtual

Definition at line 542 of file Pileup.h.

{
    myAnalyzeFuncPtr(element);
}

◆ analyzeHead()

template<class PILEUP_TYPE, class FUNC_CLASS>

void Pileup< PILEUP_TYPE, FUNC_CLASS >::analyzeHead ( )

protectedvirtual

Definition at line 549 of file Pileup.h.

{
    myAnalyzeFuncPtr(myElements[pileupHead - pileupStart]);
}

◆ flushPileup() [1/3]

template<class PILEUP_TYPE, class FUNC_CLASS>

void Pileup< PILEUP_TYPE, FUNC_CLASS >::flushPileup ( )

Done processing, flush every position that is currently being stored in the pileup.

Definition at line 368 of file Pileup.h.

{
    // while there are still entries between the head and tail, flush,
    // but no need to flush if pileupTail == -1 because in that case 
    // no entries have been added
    while ((pileupHead <= pileupTail) && (pileupTail != -1))
    {
        flushPileup(pileupHead+1);
    }
    pileupStart = pileupHead = 0;
    pileupTail = -1;
}

References flushPileup().

Referenced by ~Pileup(), flushPileup(), processAlignment(), processAlignmentRegion(), and processFile().

◆ flushPileup() [2/3]

template<class PILEUP_TYPE, class FUNC_CLASS>

void Pileup< PILEUP_TYPE, FUNC_CLASS >::flushPileup	(	int	refID,
		int	refPosition )

protectedvirtual

Definition at line 421 of file Pileup.h.

{
    // if new chromosome, flush the entire pileup.
    if(refID != myCurrentRefID)
    {
        // New chromosome, flush everything.
        flushPileup();
        myCurrentRefID = refID;
    }
    else
    {
        // on the same chromosome, so flush just up to this new position.
        flushPileup(position);
    }
}

◆ flushPileup() [3/3]

template<class PILEUP_TYPE, class FUNC_CLASS>

void Pileup< PILEUP_TYPE, FUNC_CLASS >::flushPileup ( int refPosition )

protected

Definition at line 439 of file Pileup.h.

{
    // Flush up to this new position, but no reason to flush if
    // pileupHead has not been set.
    while((pileupHead < position) && (pileupHead <= pileupTail))
    {
        analyzeHead();
 
        pileupHead++;
        
        if(pileupHead - pileupStart >= pileupWindow)
            pileupStart += pileupWindow;
    }
 
    if(pileupHead > pileupTail)
    {
        // All positions have been flushed, so reset pileup info
        pileupHead = pileupStart = 0;
        pileupTail = -1;
    }
}

◆ pileupPosition()

template<class PILEUP_TYPE, class FUNC_CLASS>

int Pileup< PILEUP_TYPE, FUNC_CLASS >::pileupPosition ( int refPosition )

protected

Definition at line 466 of file Pileup.h.

{
    // Check to see if this is the first position (pileupTail == -1)
    if(pileupTail == -1)
    {
        pileupStart = pileupHead = position;
        // This is the first time this position is being used, so
        // reset the element.
        resetElement(myElements[0], position);
        pileupTail = position;
        return(0);
    }
 
 
    if((position < pileupHead) || (position > (pileupHead + pileupWindow)))
    {
        String errorMessage =
            "Overflow on the pileup buffer: specifiedPosition = ";
        errorMessage += position;
        errorMessage += ", pileup buffer start position: ";
        errorMessage += pileupHead;
        errorMessage += ", pileup buffer end position: ";
        errorMessage += pileupHead + pileupWindow;
 
        throw std::runtime_error(errorMessage.c_str());
    }    
 
    //   int offset = position - pileupStart;
    int offset = position - pileupStart;
    
    if(offset >= pileupWindow)
    {
        offset -= pileupWindow;
    }
 
    // Check to see if position is past the end of the currently
    // setup pileup positions.
    while(position > pileupTail)
    {
        // Increment pileupTail to the next position since the current
        // pileupTail is already in use.
        ++pileupTail;
 
        // Figure out the offset for this next position.
        offset = pileupTail - pileupStart;
        if(offset >= pileupWindow)
        {
            offset -= pileupWindow;
        }
 
        // This is the first time this position is being used, so
        // reset the element.
        resetElement(myElements[offset], pileupTail);
    }
 
    return(offset);
}

◆ processAlignment()

template<class PILEUP_TYPE, class FUNC_CLASS>

void Pileup< PILEUP_TYPE, FUNC_CLASS >::processAlignment ( SamRecord & record )

virtual

Add an alignment to the pileup.

Definition at line 296 of file Pileup.h.

{
    int refPosition = record.get0BasedPosition();
    int refID = record.getReferenceID();
 
    // Flush any elements from the pileup that are prior to this record
    // since the file is sorted, we are done with those positions.
    flushPileup(refID, refPosition);
    
    // Loop through for each reference position covered by the record.
    // It is up to the PILEUP_TYPE to handle insertions/deletions, etc
    // that are related with the given reference position.
    for(; refPosition <= record.get0BasedAlignmentEnd(); ++refPosition)
    {
        addAlignmentPosition(refPosition, record);
    }
}

References flushPileup(), SamRecord::get0BasedAlignmentEnd(), SamRecord::get0BasedPosition(), and SamRecord::getReferenceID().

Referenced by processFile().

◆ processAlignmentRegion()

template<class PILEUP_TYPE, class FUNC_CLASS>

void Pileup< PILEUP_TYPE, FUNC_CLASS >::processAlignmentRegion	(	SamRecord &	record,
		int	startPos,
		int	endPos,
		PosList *	excludeList = NULL )

virtual

Add only positions that fall within the specified region of the alignment to the pileup and outside of the specified excluded positions.

Parameters

record	alignment to be added to the pileup.
startPos	0-based start position of the bases that should be added to the pileup.
endPos	0-based end position of the bases that should be added to the pileup (this position is not added). Set to -1 if there is no end position to the region.
excludeList	list of refID/positions to exclude from processing.

Definition at line 316 of file Pileup.h.

{
    int refPosition = record.get0BasedPosition();
    int refID = record.getReferenceID();
 
    // Flush any elements from the pileup that are prior to this record
    // since the file is sorted, we are done with those positions.
    flushPileup(refID, refPosition);
    
    // Check if the region starts after this reference starts.  If so,
    // we only want to start adding at the region start position.
    if(startPos > refPosition)
    {
        refPosition = startPos;
    }
 
    // Loop through for each reference position covered by the record.
    // It is up to the PILEUP_TYPE to handle insertions/deletions, etc
    // that are related with the given reference position.
    for(; refPosition <= record.get0BasedAlignmentEnd(); ++refPosition)
    {
        // Check to see if we have gone past the end of the region, in which
        // case we can stop processing this record.  Check >= since the
        // end position is not in the region.
        if((endPos != -1) && (refPosition >= endPos))
        {
            break;
        }
 
        // Check to see if this position is in the exclude list.
        bool addPos = true;
        if(excludeList != NULL)
        {
            // There is an exclude list, so lookup the position.
            if(excludeList->hasPosition(refID, refPosition))
            {
                // This position is in the exclude list, so don't add it.
                addPos = false;
            }
        }
        if(addPos)
        {
            addAlignmentPosition(refPosition, record);
        }
    }
}

References flushPileup(), SamRecord::get0BasedAlignmentEnd(), SamRecord::get0BasedPosition(), SamRecord::getReferenceID(), and PosList::hasPosition().

◆ processFile()

template<class PILEUP_TYPE, class FUNC_CLASS>

int Pileup< PILEUP_TYPE, FUNC_CLASS >::processFile	(	const std::string &	fileName,
		uint16_t	excludeFlag = 0x0704,
		uint16_t	includeFlag = 0 )

virtual

Performs a pileup on the specified file.

Parameters

excludeFlag	if specified, if any bit set in the exclude flag is set in the record's flag, it will be dropped. Defaulted to exclude: unmapped, not primary alignment failed platform/vendor quality check PCR or optical duplicate
includeFlag	if specified, every bit must be set in the record's flag for it to be included - defaulted to 0, no bits are required to be set.

Returns: 0 for success and non-zero for failure.

Definition at line 235 of file Pileup.h.

{
    SamFile samIn;
    SamFileHeader header;
    SamRecord record;
    
    if(myRefPtr != NULL)
    {
        samIn.SetReference(myRefPtr);
    }
 
    if(!samIn.OpenForRead(fileName.c_str()))
    {
        fprintf(stderr, "%s\n", samIn.GetStatusMessage());
        return(samIn.GetStatus());
    }
    
    if(!samIn.ReadHeader(header))
    {
        fprintf(stderr, "%s\n", samIn.GetStatusMessage());
        return(samIn.GetStatus());
    }
 
    // The file needs to be sorted by coordinate.
    samIn.setSortedValidation(SamFile::COORDINATE);
 
    // Iterate over all records
    while (samIn.ReadRecord(header, record))
    {
        uint16_t flag = record.getFlag();
        if(flag & excludeFlag)
        {
            // This record has an excluded flag set, 
            // so continue to the next one.
            continue;
        }
        if((flag & includeFlag) != includeFlag)
        {
            // This record does not have all required flags set, 
            // so continue to the next one.
            continue;
        }
        processAlignment(record);
    }
 
    flushPileup();
 
    int returnValue = 0;
    if(samIn.GetStatus() != SamStatus::NO_MORE_RECS)
    {
        // Failed to read a record.
        fprintf(stderr, "%s\n", samIn.GetStatusMessage());
        returnValue = samIn.GetStatus();
    }
    return(returnValue);  
}

References SamFile::COORDINATE, flushPileup(), SamRecord::getFlag(), SamFile::GetStatus(), SamFile::GetStatusMessage(), StatGenStatus::NO_MORE_RECS, SamFile::OpenForRead(), processAlignment(), SamFile::ReadHeader(), SamFile::ReadRecord(), SamFile::SetReference(), and SamFile::setSortedValidation().

◆ resetElement()

template<class PILEUP_TYPE, class FUNC_CLASS>

void Pileup< PILEUP_TYPE, FUNC_CLASS >::resetElement	(	PILEUP_TYPE &	element,
		int	position )

protectedvirtual

Definition at line 526 of file Pileup.h.

{
    element.reset(position);
}

Member Data Documentation

◆ myAnalyzeFuncPtr

template<class PILEUP_TYPE, class FUNC_CLASS = defaultPileup<PILEUP_TYPE>>

FUNC_CLASS Pileup< PILEUP_TYPE, FUNC_CLASS >::myAnalyzeFuncPtr

protected

Definition at line 119 of file Pileup.h.

◆ myCurrentRefID

template<class PILEUP_TYPE, class FUNC_CLASS = defaultPileup<PILEUP_TYPE>>

int Pileup< PILEUP_TYPE, FUNC_CLASS >::myCurrentRefID

protected

Definition at line 145 of file Pileup.h.

◆ myElements

template<class PILEUP_TYPE, class FUNC_CLASS = defaultPileup<PILEUP_TYPE>>

std::vector<PILEUP_TYPE> Pileup< PILEUP_TYPE, FUNC_CLASS >::myElements

protected

Definition at line 138 of file Pileup.h.

◆ myRefPtr

template<class PILEUP_TYPE, class FUNC_CLASS = defaultPileup<PILEUP_TYPE>>

GenomeSequence* Pileup< PILEUP_TYPE, FUNC_CLASS >::myRefPtr

protected

Definition at line 147 of file Pileup.h.

◆ pileupHead

template<class PILEUP_TYPE, class FUNC_CLASS = defaultPileup<PILEUP_TYPE>>

int Pileup< PILEUP_TYPE, FUNC_CLASS >::pileupHead

protected

Definition at line 141 of file Pileup.h.

◆ pileupStart

template<class PILEUP_TYPE, class FUNC_CLASS = defaultPileup<PILEUP_TYPE>>

int Pileup< PILEUP_TYPE, FUNC_CLASS >::pileupStart

protected

Definition at line 140 of file Pileup.h.

◆ pileupTail

template<class PILEUP_TYPE, class FUNC_CLASS = defaultPileup<PILEUP_TYPE>>

int Pileup< PILEUP_TYPE, FUNC_CLASS >::pileupTail

protected

Definition at line 142 of file Pileup.h.

◆ pileupWindow

template<class PILEUP_TYPE, class FUNC_CLASS = defaultPileup<PILEUP_TYPE>>

int Pileup< PILEUP_TYPE, FUNC_CLASS >::pileupWindow

protected

Definition at line 143 of file Pileup.h.

The documentation for this class was generated from the following file:

bam/Pileup.h

Public Member Functions

Protected Member Functions

Protected Attributes

Detailed Description

Constructor & Destructor Documentation

◆ Pileup() [1/4]

◆ Pileup() [2/4]

◆ Pileup() [3/4]

◆ Pileup() [4/4]

◆ ~Pileup()

Member Function Documentation

◆ addAlignmentPosition()

◆ addElement()

◆ analyzeElement()

◆ analyzeHead()

◆ flushPileup() [1/3]

◆ flushPileup() [2/3]

◆ flushPileup() [3/3]

◆ pileupPosition()

◆ processAlignment()

◆ processAlignmentRegion()

◆ processFile()

◆ resetElement()

Member Data Documentation

◆ myAnalyzeFuncPtr

◆ myCurrentRefID

◆ myElements

◆ myRefPtr

◆ pileupHead

◆ pileupStart

◆ pileupTail

◆ pileupWindow