SamQuerySeqWithRefHelper.cpp
00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include <stdexcept>
00019
00020 #include "SamQuerySeqWithRefHelper.h"
00021 #include "BaseUtilities.h"
00022 #include "SamFlag.h"
00023
00024 SamQuerySeqWithRefIter::SamQuerySeqWithRefIter(SamRecord& record,
00025 GenomeSequence& refSequence,
00026 bool forward)
00027 : myRecord(record),
00028 myRefSequence(refSequence),
00029 myCigar(NULL),
00030 myStartOfReadOnRefIndex(0),
00031 myQueryIndex(0),
00032 myForward(forward)
00033 {
00034 myCigar = myRecord.getCigarInfo();
00035 myStartOfReadOnRefIndex =
00036 myRefSequence.getGenomePosition(myRecord.getReferenceName()) +
00037 myRecord.get0BasedPosition();
00038
00039 if(!forward)
00040 {
00041 myQueryIndex = myRecord.getReadLength() - 1;
00042 }
00043 }
00044
00045
00046 SamQuerySeqWithRefIter::~SamQuerySeqWithRefIter()
00047 {
00048 }
00049
00050
00051
00052 bool SamQuerySeqWithRefIter::reset(bool forward)
00053 {
00054 myCigar = myRecord.getCigarInfo();
00055 if(myCigar == NULL)
00056 {
00057
00058 return(false);
00059 }
00060
00061
00062
00063 myStartOfReadOnRefIndex =
00064 myRefSequence.getGenomePosition(myRecord.getReferenceName()) +
00065 myRecord.get0BasedPosition();
00066
00067 forward = myForward;
00068
00069 if(myForward)
00070 {
00071 myQueryIndex = 0;
00072 }
00073 else
00074 {
00075
00076 myQueryIndex = myRecord.getReadLength() - 1;
00077 }
00078 return(true);
00079 }
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091 bool SamQuerySeqWithRefIter::getNextMatchMismatch(SamSingleBaseMatchInfo& matchMismatchInfo)
00092 {
00093
00094
00095 if(!SamFlag::isMapped(myRecord.getFlag()))
00096 {
00097
00098 return(false);
00099 }
00100
00101
00102 if(myCigar == NULL)
00103 {
00104
00105 throw(std::runtime_error("Cannot determine matches/mismatches since failed to retrieve the cigar"));
00106 return(false);
00107 }
00108
00109
00110
00111
00112
00113 while((myQueryIndex < myRecord.getReadLength()) &&
00114 (myQueryIndex >= 0))
00115 {
00116
00117
00118
00119 int32_t refOffset = myCigar->getRefOffset(myQueryIndex);
00120 if(refOffset == Cigar::INDEX_NA)
00121 {
00122
00123
00124
00125 nextIndex();
00126 continue;
00127 }
00128
00129
00130 char readBase = myRecord.getSequence(myQueryIndex);
00131 char refBase = myRefSequence[myStartOfReadOnRefIndex + refOffset];
00132
00133
00134
00135 if(BaseUtilities::isAmbiguous(readBase) ||
00136 BaseUtilities::isAmbiguous(refBase))
00137 {
00138
00139
00140 nextIndex();
00141 continue;
00142 }
00143
00144
00145
00146 matchMismatchInfo.setQueryIndex(myQueryIndex);
00147
00148
00149 if(BaseUtilities::areEqual(readBase, refBase))
00150 {
00151
00152 matchMismatchInfo.setType(SamSingleBaseMatchInfo::MATCH);
00153
00154 nextIndex();
00155 return(true);
00156 }
00157 else
00158 {
00159
00160 matchMismatchInfo.setType(SamSingleBaseMatchInfo::MISMATCH);
00161
00162 nextIndex();
00163 return(true);
00164 }
00165 }
00166
00167
00168 return(false);
00169 }
00170
00171
00172 void SamQuerySeqWithRefIter::nextIndex()
00173 {
00174 if(myForward)
00175 {
00176 ++myQueryIndex;
00177 }
00178 else
00179 {
00180 --myQueryIndex;
00181 }
00182 }
00183
00184
00185 SamSingleBaseMatchInfo::SamSingleBaseMatchInfo()
00186 : myType(UNKNOWN),
00187 myQueryIndex(0)
00188 {
00189 }
00190
00191
00192 SamSingleBaseMatchInfo::~SamSingleBaseMatchInfo()
00193 {
00194 }
00195
00196
00197 SamSingleBaseMatchInfo::Type SamSingleBaseMatchInfo::getType()
00198 {
00199 return(myType);
00200 }
00201
00202
00203 int32_t SamSingleBaseMatchInfo::getQueryIndex()
00204 {
00205 return(myQueryIndex);
00206 }
00207
00208
00209 void SamSingleBaseMatchInfo::setType(Type newType)
00210 {
00211 myType = newType;
00212 }
00213
00214
00215 void SamSingleBaseMatchInfo::setQueryIndex(int32_t queryIndex)
00216 {
00217 myQueryIndex = queryIndex;
00218 }