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. | |
| bool | getNextMatchMismatch (SamSingleBaseMatchInfo &matchMismatchInfo) |
| Returns information for the next position where the query and the reference match or mismatch. | |
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.
| 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
| matchMismatchInfo | return parameter with the information about the matching/mismatching base. |
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.
| 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. |
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 }
1.6.3