SamQuerySeqWithRefIter Class Reference

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

#include <SamQuerySeqWithRefHelper.h>

Collaboration diagram for SamQuerySeqWithRefIter:
Collaboration graph
[legend]

List of all members.

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.
bool getNextMatchMismatch (SamSingleBaseMatchInfo &matchMismatchInfo)
 Returns information for the next position where the query and the reference match or mismatch.

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.


Member Function Documentation

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 90 of file SamQuerySeqWithRefHelper.cpp.

References SamRecord::getFlag(), SamRecord::getReadLength(), SamRecord::getSequence(), SamRecord::NONE, SamSingleBaseMatchInfo::setQueryIndex(), and SamSingleBaseMatchInfo::setType().

00091 {
00092     // Check whether or not this read is mapped. 
00093     // If the read is not mapped, return no matches.
00094     if(!SamFlag::isMapped(myRecord.getFlag()))
00095     {
00096         // Not mapped.
00097         return(false);
00098     }
00099 
00100     // Check that the Cigar is set.
00101     if(myCigar == NULL)
00102     {
00103         // Error.
00104         throw(std::runtime_error("Cannot determine matches/mismatches since failed to retrieve the cigar"));
00105         return(false);
00106     }
00107 
00108     // Repull the read length from the record to check just in case the
00109     // record has changed length.
00110     // Loop until a match or mismatch is found as long as query index
00111     // is still within the read  (Loop is broken by a return).
00112     while((myQueryIndex < myRecord.getReadLength()) &&
00113           (myQueryIndex >= 0))
00114     {
00115         // Still more bases, look for a match/mismatch.
00116 
00117         // Get the reference offset for this read position.
00118         int32_t refOffset = myCigar->getRefOffset(myQueryIndex);
00119         if(refOffset == Cigar::INDEX_NA)
00120         {
00121             // This is either a softclip or an insertion
00122             // which do not count as a match or a mismatch, so
00123             // go to the next index.
00124             nextIndex();
00125             continue;
00126         }
00127         
00128         // Both the reference and the read have a base, so get the bases.
00129         char readBase = myRecord.getSequence(myQueryIndex, SamRecord::NONE);
00130         char refBase = myRefSequence[myStartOfReadOnRefIndex + refOffset];
00131        
00132         // If either the read or the reference bases are unknown, then
00133         // it does not count as a match or a mismatch.
00134         if(BaseUtilities::isAmbiguous(readBase) || 
00135            BaseUtilities::isAmbiguous(refBase))
00136         {
00137             // Either the reference base or the read base are unknown,
00138             // so skip this position.
00139             nextIndex();
00140             continue;
00141         }
00142 
00143         // Both the read & the reference have a known base, so it is either
00144         // a match or a mismatch.
00145         matchMismatchInfo.setQueryIndex(myQueryIndex);
00146 
00147         // Check if they are equal.
00148         if(BaseUtilities::areEqual(readBase, refBase))
00149         {
00150             // Match.
00151             matchMismatchInfo.setType(SamSingleBaseMatchInfo::MATCH);
00152             // Increment the query index to the next position.
00153             nextIndex();
00154             return(true);
00155         }
00156         else
00157         {
00158             // Mismatch
00159             matchMismatchInfo.setType(SamSingleBaseMatchInfo::MISMATCH);
00160             // Increment the query index to the next position.
00161             nextIndex();
00162             return(true);
00163         }
00164     }
00165 
00166     // No matches or mismatches were found, so return false.
00167     return(false);
00168 }

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 52 of file SamQuerySeqWithRefHelper.cpp.

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

00053 {
00054     myCigar = myRecord.getCigarInfo();
00055     if(myCigar == NULL)
00056     {
00057         // Failed to get Cigar.
00058         return(false);
00059     }
00060 
00061     // Get where the position of where this read starts as mapped to the 
00062     // reference.
00063     myStartOfReadOnRefIndex = 
00064         myRefSequence.getGenomePosition(myRecord.getReferenceName()) +
00065         myRecord.get0BasedPosition();
00066 
00067     myForward = forward;
00068 
00069     if(myForward)
00070     {
00071         myQueryIndex = 0;
00072     }
00073     else
00074     {
00075         // reverse, so start at the last entry.
00076         myQueryIndex = myRecord.getReadLength() - 1;
00077     }
00078     return(true);
00079 }


The documentation for this class was generated from the following files:
Generated on Tue Mar 22 22:50:27 2011 for StatGen Software by  doxygen 1.6.3