Skip to content

Commit

Permalink
incorporate changes for co-linear consensus mechanism
Browse files Browse the repository at this point in the history
  • Loading branch information
@rob-p
rob-p committed Apr 16, 2016
1 parent bc7d801 commit 617e48b
Show file tree
Hide file tree
Showing 5 changed files with 194 additions and 131 deletions.
12 changes: 7 additions & 5 deletions include/HitManager.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,9 +72,10 @@ namespace rapmap {

template <typename RapMapIndexT>
void intersectSAIntervalWithOutput(SAIntervalHit<typename RapMapIndexT::IndexType>& h,
RapMapIndexT& rmi,
uint32_t intervalCounter,
SAHitMap& outHits);
RapMapIndexT& rmi,
uint32_t intervalCounter,
SAHitMap& outHits);


template <typename RapMapIndexT>
void intersectSAIntervalWithOutput2(SAIntervalHit<typename RapMapIndexT::IndexType>& h,
Expand All @@ -93,8 +94,9 @@ namespace rapmap {

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

template <typename RapMapIndexT>
std::vector<ProcessedSAHit> intersectSAHits2(
Expand Down
42 changes: 41 additions & 1 deletion include/RapMapUtils.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -369,9 +369,49 @@ namespace rapmap {
tqvec.emplace_back(txpPosIn, queryPosIn, queryRCIn);
}

/**
* This enforces a more stringent consistency check on
* the hits for this transcript. The hits must be co-linear
* with respect to the query and target.
*
* input: numToCheck --- the number of hits to check in sorted order
* hits after the last of these need not be consistent.
* return: numToCheck if the first numToCheck hits are consistent;
* -1 otherwise
**/
int32_t checkConsistent(int32_t numToCheck) {
auto numHits = tqvec.size();

// special case for only 1 or two hits (common)
if (numHits == 1) {
return numToCheck;
} else if (numHits == 2) {
auto& h1 = (tqvec[0].queryPos < tqvec[1].queryPos) ? tqvec[0] : tqvec[1];
auto& h2 = (tqvec[0].queryPos < tqvec[1].queryPos) ? tqvec[1] : tqvec[2];
return (h2.pos > h1.pos) ? (numToCheck) : -1;
} else {
// first, sort by query position
std::sort(tqvec.begin(), tqvec.end(),
[](const SATxpQueryPos& q1, const SATxpQueryPos& q2) -> bool {
return q1.queryPos < q2.queryPos;
});

int32_t lastRefPos{std::numeric_limits<int32_t>::min()};
for (size_t i = 0; i < numToCheck; ++i) {
int32_t refPos = static_cast<int32_t>(tqvec[i].pos);
if (refPos > lastRefPos) {
lastRefPos = refPos;
} else {
return i;
}
}
return numToCheck;
}
}

uint32_t tid;
std::vector<SATxpQueryPos> tqvec;
bool active;
bool active;
uint32_t numActive;
};

Expand Down
17 changes: 9 additions & 8 deletions include/SACollector.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,10 +16,11 @@ class SACollector {

SACollector(RapMapIndexT* rmi) : rmi_(rmi) {}
bool operator()(std::string& read,
std::vector<rapmap::utils::QuasiAlignment>& hits,
SASearcher<RapMapIndexT>& saSearcher,
rapmap::utils::MateStatus mateStatus,
bool strictCheck=false) {
std::vector<rapmap::utils::QuasiAlignment>& hits,
SASearcher<RapMapIndexT>& saSearcher,
rapmap::utils::MateStatus mateStatus,
bool strictCheck=false,
bool consistentHits=false) {

using QuasiAlignment = rapmap::utils::QuasiAlignment;
using MateStatus = rapmap::utils::MateStatus;
Expand Down Expand Up @@ -483,10 +484,10 @@ class SACollector {
auto fwdHitsStart = hits.size();
// If we had > 1 forward hit
if (fwdSAInts.size() > 1) {
auto processedHits = rapmap::hit_manager::intersectSAHits(fwdSAInts, *rmi_);
rapmap::hit_manager::collectHitsSimpleSA(processedHits, readLen, maxDist, hits, mateStatus);
auto processedHits = rapmap::hit_manager::intersectSAHits(fwdSAInts, *rmi_, consistentHits);
rapmap::hit_manager::collectHitsSimpleSA(processedHits, readLen, maxDist, hits, mateStatus);
} else if (fwdSAInts.size() == 1) { // only 1 hit!
auto& saIntervalHit = fwdSAInts.front();
auto& saIntervalHit = fwdSAInts.front();
auto initialSize = hits.size();
for (OffsetT i = saIntervalHit.begin; i != saIntervalHit.end; ++i) {
auto globalPos = SA[i];
Expand Down Expand Up @@ -520,7 +521,7 @@ class SACollector {
auto rcHitsStart = fwdHitsEnd;
// If we had > 1 rc hit
if (rcSAInts.size() > 1) {
auto processedHits = rapmap::hit_manager::intersectSAHits(rcSAInts, *rmi_);
auto processedHits = rapmap::hit_manager::intersectSAHits(rcSAInts, *rmi_, consistentHits);
rapmap::hit_manager::collectHitsSimpleSA(processedHits, readLen, maxDist, hits, mateStatus);
} else if (rcSAInts.size() == 1) { // only 1 hit!
auto& saIntervalHit = rcSAInts.front();
Expand Down
48 changes: 28 additions & 20 deletions src/HitManager.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -346,9 +346,9 @@ namespace rapmap {

template <typename RapMapIndexT>
void intersectSAIntervalWithOutput(SAIntervalHit<typename RapMapIndexT::IndexType>& h,
RapMapIndexT& rmi,
uint32_t intervalCounter,
SAHitMap& outHits) {
RapMapIndexT& rmi,
uint32_t intervalCounter,
SAHitMap& outHits) {
using OffsetT = typename RapMapIndexT::IndexType;
// Convenient bindings for variables we'll use
auto& SA = rmi.SA;
Expand Down Expand Up @@ -577,7 +577,8 @@ namespace rapmap {
template <typename RapMapIndexT>
SAHitMap intersectSAHits(
std::vector<SAIntervalHit<typename RapMapIndexT::IndexType>>& inHits,
RapMapIndexT& rmi
RapMapIndexT& rmi,
bool strictFilter
) {
using OffsetT = typename RapMapIndexT::IndexType;
// Each inHit is a SAIntervalHit structure that contains
Expand All @@ -592,7 +593,7 @@ namespace rapmap {
// with less than 2 hits.
SAHitMap outHits;
if (inHits.size() < 2) {
std::cerr << "intersectHitsSA() called with < 2 k-mer "
std::cerr << "intersectHitsSA() called with < 2 hits "
" hits; this shouldn't happen\n";
return outHits;
}
Expand Down Expand Up @@ -626,20 +627,21 @@ namespace rapmap {

// Now intersect everything in inHits (apart from minHits)
// to get the final set of mapping info.
size_t intervalCounter{2};
size_t intervalCounter{2};
for (auto& h : inHits) {
if (&h != minHit) { // don't intersect minHit with itself
intersectSAIntervalWithOutput(h, rmi, intervalCounter, outHits);
++intervalCounter;
++intervalCounter;
}
}

size_t requiredNumHits = inHits.size();
// Mark as active any transcripts with the required number of hits.
for (auto it = outHits.begin(); it != outHits.end(); ++it) {
if (it->second.numActive >= requiredNumHits) {
it->second.active = true;
}
bool enoughHits = (it->second.numActive >= requiredNumHits);
it->second.active = (strictFilter) ?
(enoughHits and it->second.checkConsistent(requiredNumHits)) :
(enoughHits);
}
return outHits;
}
Expand All @@ -657,36 +659,42 @@ namespace rapmap {

template
void intersectSAIntervalWithOutput<SAIndex32BitDense>(SAIntervalHit<int32_t>& h,
SAIndex32BitDense& rmi, uint32_t intervalCounter, SAHitMap& outHits);
SAIndex32BitDense& rmi,
uint32_t intervalCounter,
SAHitMap& outHits);

template
void intersectSAIntervalWithOutput<SAIndex64BitDense>(SAIntervalHit<int64_t>& h,
SAIndex64BitDense& rmi, uint32_t intervalCounter, SAHitMap& outHits);
SAIndex64BitDense& rmi,
uint32_t intervalCounter,
SAHitMap& outHits);

template
SAHitMap intersectSAHits<SAIndex32BitDense>(std::vector<SAIntervalHit<int32_t>>& inHits,
SAIndex32BitDense& rmi);
SAIndex32BitDense& rmi, bool strictFilter);

template
SAHitMap intersectSAHits<SAIndex64BitDense>(std::vector<SAIntervalHit<int64_t>>& inHits,
SAIndex64BitDense& rmi);
SAIndex64BitDense& rmi, bool strictFilter);

template
void intersectSAIntervalWithOutput<SAIndex32BitPerfect>(SAIntervalHit<int32_t>& h,
SAIndex32BitPerfect& rmi, uint32_t intervalCounter, SAHitMap& outHits);
SAIndex32BitPerfect& rmi,
uint32_t intervalCounter,
SAHitMap& outHits);

template
void intersectSAIntervalWithOutput<SAIndex64BitPerfect>(SAIntervalHit<int64_t>& h,
SAIndex64BitPerfect& rmi, uint32_t intervalCounter, SAHitMap& outHits);
SAIndex64BitPerfect& rmi,
uint32_t intervalCounter,
SAHitMap& outHits);

template
SAHitMap intersectSAHits<SAIndex32BitPerfect>(std::vector<SAIntervalHit<int32_t>>& inHits,
SAIndex32BitPerfect& rmi);
SAIndex32BitPerfect& rmi, bool strictFilter);

template
SAHitMap intersectSAHits<SAIndex64BitPerfect>(std::vector<SAIntervalHit<int64_t>>& inHits,
SAIndex64BitPerfect& rmi);


SAIndex64BitPerfect& rmi, bool strictFilter);
}
}
Loading

0 comments on commit 617e48b

Please sign in to comment.