PileupElementBaseQual.cpp
/*
* Copyright (C) 2010 Regents of the University of Michigan
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <stdexcept>
#include "PileupElementBaseQual.h"
PileupElementBaseQual::PileupElementBaseQual()
: PileupElement(),
myBases(NULL),
myQualities(NULL),
myAllocatedSize(0),
myIndex(-1),
myAddDelAsBase(false)
{
myAllocatedSize = 1024;
myBases = (char*)malloc(myAllocatedSize + 1);
myQualities = (char*)malloc(myAllocatedSize + 1);
if((myBases == NULL ) || (myQualities == NULL))
{
// TODO, check for malloc failures.
std::cerr << "Failed Memory Allocation\n";
}
}
// NOTE that this method does not actually copy, it just resets.
PileupElementBaseQual::PileupElementBaseQual(const PileupElementBaseQual& q)
: PileupElement(),
myBases(NULL),
myQualities(NULL),
myAllocatedSize(0),
myIndex(-1)
{
myAllocatedSize = 1024;
myBases = (char*)malloc(myAllocatedSize + 1);
myQualities = (char*)malloc(myAllocatedSize + 1);
myAddDelAsBase = q.myAddDelAsBase;
if((myBases == NULL ) || (myQualities == NULL))
{
// TODO, check for malloc failures.
std::cerr << "Failed Memory Allocation\n";
}
}
PileupElementBaseQual::~PileupElementBaseQual()
{
if(myBases != NULL)
{
free(myBases);
myBases = NULL;
}
if(myQualities != NULL)
{
free(myQualities);
myQualities = NULL;
}
}
// Add an entry to this pileup element.
void PileupElementBaseQual::addEntry(SamRecord& record)
{
// Call the base class:
PileupElement::addEntry(record);
// Increment the index
++myIndex;
// if the index has gone beyond the allocated space, double the size.
if(myIndex >= myAllocatedSize)
{
char* tempBuffer = (char*)realloc(myBases, myAllocatedSize * 2);
if(tempBuffer == NULL)
{
std::cerr << "Memory Allocation Failure\n";
// TODO
return;
}
myBases = tempBuffer;
tempBuffer = (char*)realloc(myQualities, myAllocatedSize * 2);
if(tempBuffer == NULL)
{
std::cerr << "Memory Allocation Failure\n";
// TODO
return;
}
myQualities = tempBuffer;
myAllocatedSize = myAllocatedSize * 2;
}
Cigar* cigar = record.getCigarInfo();
if(cigar == NULL)
{
throw std::runtime_error("Failed to retrieve cigar info from the record.");
}
int32_t readIndex =
cigar->getQueryIndex(getRefPosition(), record.get0BasedPosition());
// If the readPosition is N/A, this is a deletion.
if(readIndex != CigarRoller::INDEX_NA)
{
char base = record.getSequence(readIndex);
char qual = record.getQuality(readIndex);
if(qual == UNSET_QUAL)
{
qual = ' ';
}
myBases[myIndex] = base;
myQualities[myIndex] = qual;
}
else if(myAddDelAsBase)
{
// This version adds deletions as bases.
myBases[myIndex] = '-';
myQualities[myIndex] = '0';
}
else
{
// Do not add a deletion.
// Did not add any entries, so decrement the index counter since the
// index was not used.
--myIndex;
}
}
void PileupElementBaseQual::analyze()
{
if(getRefPosition() != UNSET_POSITION)
{
myBases[myIndex+1] = '\0';
myQualities[myIndex+1] = '\0';
std::cout << getChromosome() << "\t" << getRefPosition() << "\tN\t" << myIndex+1 << "\t";
std::cout << myBases << "\t";
std::cout << myQualities;
std::cout << "\n";
}
}
void PileupElementBaseQual::reset(int refPosition)
{
// Call the base class.
PileupElement::reset(refPosition);
myIndex = -1;
}
Cigar
This class represents the CIGAR without any methods to set the cigar (see CigarRoller for that).
Definition: Cigar.h:84
Cigar::INDEX_NA
static const int32_t INDEX_NA
Value associated with an index that is not applicable/does not exist, used for converting between que...
Definition: Cigar.h:492
Cigar::getQueryIndex
int32_t getQueryIndex(int32_t refOffset)
Return the query index associated with the specified reference offset or INDEX_NA based on this cigar...
Definition: Cigar.cpp:202
PileupElementBaseQual
This class inherits from the base class and stores base and qualities.
Definition: PileupElementBaseQual.h:26
PileupElementBaseQual::reset
virtual void reset(int32_t refPosition)
Resets the entry, setting the new position associated with this element.
Definition: PileupElementBaseQual.cpp:157
PileupElementBaseQual::analyze
virtual void analyze()
Perform the analysis associated with this class.
Definition: PileupElementBaseQual.cpp:144
PileupElementBaseQual::addEntry
virtual void addEntry(SamRecord &record)
Add an entry to this pileup element.
Definition: PileupElementBaseQual.cpp:76
PileupElement
This is a base class pileup component, representing the information for one reference position.
Definition: PileupElement.h:27
PileupElement::getRefPosition
int32_t getRefPosition() const
Get the reference position stored in this element.
Definition: PileupElement.h:54
PileupElement::getChromosome
const char * getChromosome() const
Get the chromosome name stored in this element.
Definition: PileupElement.h:51
PileupElement::addEntry
virtual void addEntry(SamRecord &record)
Add an entry to this pileup element.
Definition: PileupElement.cpp:44
PileupElement::reset
virtual void reset(int32_t refPosition)
Resets the entry, setting the new position associated with this element.
Definition: PileupElement.cpp:67
SamRecord
Class providing an easy to use interface to get/set/operate on the fields in a SAM/BAM record.
Definition: SamRecord.h:52
SamRecord::getCigarInfo
Cigar * getCigarInfo()
Returns a pointer to the Cigar object associated with this record.
Definition: SamRecord.cpp:1836
SamRecord::get0BasedPosition
int32_t get0BasedPosition()
Get the 0-based(BAM) leftmost position of the record.
Definition: SamRecord.cpp:1319
SamRecord::getQuality
const char * getQuality()
Returns the SAM formatted quality string (QUAL).
Definition: SamRecord.cpp:1638
SamRecord::getSequence
const char * getSequence()
Returns the SAM formatted sequence string (SEQ), translating the base as specified by setSequenceTran...
Definition: SamRecord.cpp:1568