SamQuerySeqWithRefIter Class Reference

Iterates through the query and compare with reference. More...

`#include <SamQuerySeqWithRefHelper.h>`

## Public Member Functions

| SamQuerySeqWithRefIter (SamRecord &record, GenomeSequence &refSequence, bool forward=true)
bool | reset (bool forward=true)
| Reset to start at the beginning of the record. More...
bool | getNextMatchMismatch (SamSingleBaseMatchInfo &matchMismatchInfo)
| Returns information for the next position where the query and the reference match or mismatch. More...

## Detailed Description

Iterates through the query and compare with reference.

NOTE: References to the GenomeSequence and SamRecord are stored, the objects are not copied, so they must remain valid as long as this class is used.

Definition at line 58 of file SamQuerySeqWithRefHelper.h.

## Constructor & Destructor Documentation

## SamQuerySeqWithRefIter()

SamQuerySeqWithRefIter::SamQuerySeqWithRefIter | ( | SamRecord & | record,
|  | GenomeSequence & | refSequence,
|  | bool | forward = `true`
| ) |  |

Definition at line 24 of file SamQuerySeqWithRefHelper.cpp.

: myRecord(record),
myRefSequence(refSequence),
myCigar(NULL),
myStartOfReadOnRefIndex(INVALID_GENOME_INDEX),
myQueryIndex(0),
myForward(forward)
{
myCigar = myRecord.getCigarInfo();
myStartOfReadOnRefIndex =
refSequence.getGenomePosition(myRecord.getReferenceName());
if(myStartOfReadOnRefIndex != INVALID_GENOME_INDEX)
{
// This reference name was found in the reference file, so
// add the start position.
myStartOfReadOnRefIndex += myRecord.get0BasedPosition();
}
if(!forward)
{
myQueryIndex = myRecord.getReadLength() - 1;
}
}
GenomeSequence::getGenomePosition
genomeIndex_t getGenomePosition(const char *chromosomeName, unsigned int chromosomeIndex) const
given a chromosome name and position, return the genome position
Definition: GenomeSequence.cpp:779
SamRecord::getReferenceName
const char * getReferenceName()
Get the reference sequence name (RNAME) of the record.
Definition: SamRecord.cpp:1298
SamRecord::getCigarInfo
Cigar * getCigarInfo()
Returns a pointer to the Cigar object associated with this record.
Definition: SamRecord.cpp:1836
SamRecord::getReadLength
int32_t getReadLength()
Get the length of the read.
Definition: SamRecord.cpp:1391
SamRecord::get0BasedPosition
int32_t get0BasedPosition()
Get the 0-based(BAM) leftmost position of the record.
Definition: SamRecord.cpp:1319

## ~SamQuerySeqWithRefIter()

SamQuerySeqWithRefIter::~SamQuerySeqWithRefIter | ( |  | ) |

virtual

Definition at line 52 of file SamQuerySeqWithRefHelper.cpp.

{
}

## Member Function Documentation

## getNextMatchMismatch()

bool SamQuerySeqWithRefIter::getNextMatchMismatch | ( | SamSingleBaseMatchInfo & | matchMismatchInfo | ) |

Returns information for the next position where the query and the reference match or mismatch.

To be a match or mismatch, both the query and reference must have a base that is not 'N'. This means: insertions and deletions are not mismatches or matches. 'N' bases are not matches or mismatches

Parameters

  matchMismatchInfo | return parameter with the information about the matching/mismatching base.

Returns
  true if there was another match/mismatch (matchMismatchInfo was set); false if not.

Definition at line 102 of file SamQuerySeqWithRefHelper.cpp.

{
// Check whether or not this read is mapped.
// If the read is not mapped, return no matches.
if(!SamFlag::isMapped(myRecord.getFlag()))
{
// Not mapped.
return(false);
}
// Check that the Cigar is set.
if(myCigar == NULL)
{
// Error.
throw(std::runtime_error("Cannot determine matches/mismatches since failed to retrieve the cigar"));
return(false);
}
// If myStartOfReadOnRefIndex is the default (unset) value, then
// the reference was not found, so return false, no matches/mismatches.
if(myStartOfReadOnRefIndex == INVALID_GENOME_INDEX)
{
// This reference name was not found in the reference file, so just
// return no matches/mismatches.
return(false);
}
// Repull the read length from the record to check just in case the
// record has changed length.
// Loop until a match or mismatch is found as long as query index
// is still within the read (Loop is broken by a return).
while((myQueryIndex < myRecord.getReadLength()) &&
(myQueryIndex >= 0))
{
// Still more bases, look for a match/mismatch.
// Get the reference offset for this read position.
int32_t refOffset = myCigar->getRefOffset(myQueryIndex);
if(refOffset == Cigar::INDEX_NA)
{
// This is either a softclip or an insertion
// which do not count as a match or a mismatch, so
// go to the next index.
nextIndex();
continue;
}
// Both the reference and the read have a base, so get the bases.
char readBase = myRecord.getSequence(myQueryIndex, SamRecord::NONE);
char refBase = myRefSequence[myStartOfReadOnRefIndex + refOffset];
// If either the read or the reference bases are unknown, then
// it does not count as a match or a mismatch.
if(BaseUtilities::isAmbiguous(readBase) ||
BaseUtilities::isAmbiguous(refBase))
{
// Either the reference base or the read base are unknown,
// so skip this position.
nextIndex();
continue;
}
// Both the read & the reference have a known base, so it is either
// a match or a mismatch.
matchMismatchInfo.setQueryIndex(myQueryIndex);
// Check if they are equal.
if(BaseUtilities::areEqual(readBase, refBase))
{
// Match.
matchMismatchInfo.setType(SamSingleBaseMatchInfo::MATCH);
// Increment the query index to the next position.
nextIndex();
return(true);
}
else
{
// Mismatch
matchMismatchInfo.setType(SamSingleBaseMatchInfo::MISMATCH);
// Increment the query index to the next position.
nextIndex();
return(true);
}
}
// No matches or mismatches were found, so return false.
return(false);
}
BaseUtilities::isAmbiguous
static bool isAmbiguous(char base)
Returns whether or not the specified bases is an indicator for ambiguity.
Definition: BaseUtilities.cpp:23
BaseUtilities::areEqual
static bool areEqual(char base1, char base2)
Returns whether or not two bases are equal (case insensitive), if one of the bases is '=',...
Definition: BaseUtilities.cpp:39
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::getRefOffset
int32_t getRefOffset(int32_t queryIndex)
Return the reference offset associated with the specified query index or INDEX_NA based on this cigar...
Definition: Cigar.cpp:187
SamRecord::NONE
@ NONE
Leave the sequence as is.
Definition: SamRecord.h:58
SamRecord::getFlag
uint16_t getFlag()
Get the flag (FLAG).
Definition: SamRecord.cpp:1384
SamRecord::getSequence
const char * getSequence()
Returns the SAM formatted sequence string (SEQ), translating the base as specified by setSequenceTran...
Definition: SamRecord.cpp:1568
SamSingleBaseMatchInfo::setQueryIndex
void setQueryIndex(int32_t queryIndex)
Set the query index for this object.
Definition: SamQuerySeqWithRefHelper.cpp:236
SamSingleBaseMatchInfo::setType
void setType(Type newType)
Set the type (match/mismatch/unkown) for this object.
Definition: SamQuerySeqWithRefHelper.cpp:230

References BaseUtilities::areEqual(), SamRecord::getFlag(), SamRecord::getReadLength(), Cigar::getRefOffset(), SamRecord::getSequence(), Cigar::INDEX_NA, BaseUtilities::isAmbiguous(), SamRecord::NONE, SamSingleBaseMatchInfo::setQueryIndex(), and SamSingleBaseMatchInfo::setType().

Referenced by SamFilter::clipOnMismatchThreshold(), and SamFilter::sumMismatchQuality().

## reset()

bool SamQuerySeqWithRefIter::reset | ( | bool | forward = `true` | ) |

Reset to start at the beginning of the record.

This will re-read values from SamRecord, so can be used if it has changed to contain information for a new record.

Parameters

  forward | true means to start from the beginning and go to the end; false means to start from the end and go to the beginning.

Returns
  true if successfully reset; false if failed to read the Cigar.

Definition at line 58 of file SamQuerySeqWithRefHelper.cpp.

{
myCigar = myRecord.getCigarInfo();
if(myCigar == NULL)
{
// Failed to get Cigar.
return(false);
}
// Get where the position of where this read starts as mapped to the
// reference.
myStartOfReadOnRefIndex =
myRefSequence.getGenomePosition(myRecord.getReferenceName());
if(myStartOfReadOnRefIndex != INVALID_GENOME_INDEX)
{
// This reference name was found in the reference file, so
// add the start position.
myStartOfReadOnRefIndex += myRecord.get0BasedPosition();
}
myForward = forward;
if(myForward)
{
myQueryIndex = 0;
}
else
{
// reverse, so start at the last entry.
myQueryIndex = myRecord.getReadLength() - 1;
}
return(true);
}

References SamRecord::get0BasedPosition(), SamRecord::getCigarInfo(), GenomeSequence::getGenomePosition(), SamRecord::getReadLength(), and SamRecord::getReferenceName().

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

- bam/SamQuerySeqWithRefHelper.h
- bam/SamQuerySeqWithRefHelper.cpp