Skip to content

Commit

Permalink
Address issue #46.
Browse files Browse the repository at this point in the history 
In the SAM output of RapMap, there was an implicit assumption that
if a read is unmapped, then the read is written down in the forward
orientation.  That is, if the read unmapped flag is set (0x4), then
the flag 0x10 was assumed to have an invalid (random) state ---
the unmapped read is always recorded in its forward orientation.
However, while the SAM spec does dictate that many of the flags can
be ignored (i.e. are not guaranteed to be valid) if 0x4 is set, it
actually doesn't seem to make this claim about 0x10.  These changes
should assure that, for an orphaned pair, the mateFwd flag will be
set.
  • Loading branch information
Rob Patro committed Feb 1, 2019
1 parent 23ec452 commit 836f85c
Show file tree
Hide file tree
Showing 2 changed files with 9 additions and 124 deletions.
18 changes: 1 addition & 17 deletions include/HitManager.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,35 +114,19 @@ namespace rapmap {
uint32_t intervalCounter,
int32_t maxSlack,
SAHitMap& outHits);


template <typename RapMapIndexT>
void intersectSAIntervalWithOutput2(SAIntervalHit<typename RapMapIndexT::IndexType>& h,
RapMapIndexT& rmi,
SAProcessedHitVec& outStructs);

/*
void intersectSAIntervalWithOutput3(SAIntervalHit& h,
RapMapSAIndex& rmi,
SAProcessedHitVec& outHits);
*/

//std::vector<ProcessedHit> intersectHits(
// std::vector<HitInfo>& inHits,
// RapMapIndex& rmi);

template <typename RapMapIndexT>
template <typename RapMapIndexT>
SAHitMap intersectSAHits(
SAIntervalVector<SAIntervalHit<typename RapMapIndexT::IndexType>>& inHits,
RapMapIndexT& rmi,
size_t readLen,
bool strictFilter=false);

template <typename RapMapIndexT>
std::vector<ProcessedSAHit> intersectSAHits2(
std::vector<SAIntervalHit<typename RapMapIndexT::IndexType>>& inHits,
RapMapIndexT& rmi);

template <typename RapMapIndexT>
void hitsToMappingsSimple(RapMapIndexT& rmi,
rapmap::utils::MappingConfig& mc,
Expand Down
115 changes: 8 additions & 107 deletions src/HitManager.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,6 +72,7 @@ namespace rapmap {
// determine forward
hits.emplace_back(tid, hitPos, isFwd, readLen);
hits.back().mateStatus = mateStatus;
hits.back().mateIsFwd = true;
}

return true;
Expand Down Expand Up @@ -262,6 +263,7 @@ namespace rapmap {
int32_t hitPos = (*posIt)->pos - (*posIt)->queryPos;
hits.emplace_back(tid, hitPos, isFwd, readLen);
auto& currHit = hits.back();
currHit.mateIsFwd = true; // address issue #46
currHit.setChainScore(bestScore);
currHit.mateStatus = mateStatus;
currHit.allPositions.push_back(hitPos);
Expand Down Expand Up @@ -313,119 +315,16 @@ namespace rapmap {
int32_t hitPos = minPosIt->pos - minPosIt->queryPos;
bool isFwd = !hitRC;
hits.emplace_back(tid, hitPos, isFwd, readLen);
hits.back().mateStatus = mateStatus;
hits.back().allPositions.push_back(hitPos);
auto& currHit = hits.back();
currHit.mateIsFwd = true; // address issue #46
currHit.mateStatus = mateStatus;
currHit.allPositions.push_back(hitPos);
}
}
}
// if SAHitMap is sorted, no need to sort here
/*
std::sort(hits.begin() + startOffset, hits.end(),
[](const QuasiAlignment& a, const QuasiAlignment& b) -> bool {
return a.tid < b.tid;
});
*/
return true;
}


// Return hits from processedHits where position constraints
// match maxDist
bool collectHitsSimpleSA2(std::vector<ProcessedSAHit>& processedHits,
uint32_t readLen,
uint32_t maxDist,
std::vector<QuasiAlignment>& hits,
MateStatus mateStatus){
//bool foundHit{false};

// One processed hit per transcript
for (auto& ph : processedHits) {
// If this is an *active* position list
if (ph.active) {
auto tid = ph.tid;
auto minPosIt =
std::min_element(ph.tqvec.begin(),
ph.tqvec.end(),
[](const SATxpQueryPos& a, const SATxpQueryPos& b) -> bool {
return a.pos < b.pos;
});

bool hitRC = minPosIt->queryRC;
int32_t hitPos = minPosIt->pos - minPosIt->queryPos;
bool isFwd = !hitRC;
hits.emplace_back(tid, hitPos, isFwd, readLen);
hits.back().mateStatus = mateStatus;
}
}
return true;
}




// Intersects the hit h2 with outHits.
// This will modify outHits so that the tqvec field of the
// entries in outHits that are labeled by the transcripts in
// which h2 appears will have an iterator to the beginning of
// the position list for h2.
/*void intersectWithOutput(HitInfo& h2, RapMapIndex& rmi,
std::vector<ProcessedHit>& outHits) {
// Convenient bindings for variables we'll use
auto& eqClasses = rmi.eqClassList;
auto& eqClassLabels = rmi.eqLabelList;
auto& posList = rmi.posList;
// Iterator to the beginning and end of the output hits
auto outHitIt = outHits.begin();
auto outHitEnd = outHits.end();
// Equiv. class for h2
auto& eqClassRight = eqClasses[h2.kinfo->eqId];
// Iterator into, length of and end of the positon list for h2
auto rightPosIt = posList.begin() + h2.kinfo->offset;
auto rightPosLen = h2.kinfo->count;
// auto rightPosEnd = rightPosIt + rightPosLen;
// Iterator into, length of and end of the transcript list for h2
auto rightTxpIt = eqClassLabels.begin() + eqClassRight.txpListStart;
auto rightTxpListLen = eqClassRight.txpListLen;
auto rightTxpEnd = rightTxpIt + rightTxpListLen;
auto rightQueryPos = h2.queryPos;
auto rightQueryRC = h2.queryRC;
PositionListHelper rightPosHelper(rightPosIt, posList.end());
uint32_t leftTxp, rightTxp;
while (outHitIt != outHitEnd and rightTxpIt != rightTxpEnd) {
// Get the current transcript ID for the left and right eq class
leftTxp = outHitIt->tid;
rightTxp = *rightTxpIt;
// If we need to advance the left txp, do it
if (leftTxp < rightTxp) {
// Advance to the next transcript in the
// equivalence class label
++outHitIt;
} else {
// If the transcripts are equal (i.e. leftTxp >= rightTxp and !(rightTxp < leftTxp))
// Then see if there are any hits here.
if (!(rightTxp < leftTxp)) {
// Add the position list iterator and query pos for the
// hit from h2 to the back of outHits' tqvec.
outHitIt->tqvec.emplace_back(rightPosHelper, rightQueryPos, rightQueryRC);
++outHitIt;
}
// advance the hit we're intersecting to the next transcript
rightPosHelper.advanceToNextTranscript();
// Advance the right transcript id regardless of whether
// we found a hit or not.
++rightTxpIt;
}
}
}
*/

/** from http://en.cppreference.com/w/cpp/algorithm/lower_bound **/
template <typename ForwardIt>
ForwardIt binarySearch(
Expand Down Expand Up @@ -830,6 +729,8 @@ namespace rapmap {
int32_t hitPos = pos - saIntervalHit.queryPos;
outHits.emplace_back(txpID, hitPos, isFw, readLen);
auto& lastHit = outHits.back();
lastHit.mateIsFwd = true; // address issue #46

//lastHit.queryOffset = saIntervalHit.queryPos;
lastHit.mateStatus = mateStatus;
lastHit.allPositions.push_back(hitPos);
Expand Down

0 comments on commit 836f85c

Please sign in to comment.