From 1b4ecc010b00c6c8ccd25356301fe1f22aee1e69 Mon Sep 17 00:00:00 2001
From: Aaron <aaronquinlan@gmail.com>
Date: Mon, 17 Sep 2012 16:56:57 -0400
Subject: [PATCH] fix corner case in chrom sweep. addresses GC issue #141
---
src/intersectBed/intersectBed.cpp | 20 +++++++++++++-------
src/utils/chromsweep/chromsweep.cpp | 12 ++++++++----
2 files changed, 21 insertions(+), 11 deletions(-)
diff --git a/src/intersectBed/intersectBed.cpp b/src/intersectBed/intersectBed.cpp
index 8393ea1e..3382b2df 100644
--- a/src/intersectBed/intersectBed.cpp
+++ b/src/intersectBed/intersectBed.cpp
@@ -246,14 +246,17 @@ void BedIntersect::IntersectBed() {
// treat the BED as a single "block"
if (_obeySplits == false)
FindOverlaps(a, hits);
- // split the BED12 into blocks and look for overlaps in each discrete block
+ // split the BED12 into blocks and look for
+ // overlaps in each discrete block
else {
- // find the hits that overlap with the full span of the blocked BED
+ // find the hits that overlap with the
+ // full span of the blocked BED
_bedB->allHits(a.chrom, a.start, a.end, a.strand,
hits, _sameStrand, _diffStrand,
0.0, false);
// break a into discrete blocks, as we need to
- // measure overlap with the individual blocks, not the full span.
+ // measure overlap with the individual blocks,
+ // not the full span.
bedVector a_blocks;
GetBedBlocks(a, a_blocks);
// find the overlaps between the block in A and B
@@ -275,12 +278,14 @@ void BedIntersect::IntersectBed() {
pair<BED, vector<BED> > hit_set;
hit_set.second.reserve(10000);
while (sweep.Next(hit_set)) {
- if (_obeySplits == false)
- processHits(hit_set.first, hit_set.second);
+ if (_obeySplits == false) {
+ processHits(hit_set.first, hit_set.second);
+ }
else {
bedVector a_blocks;
GetBedBlocks(hit_set.first, a_blocks);
- FindBlockedOverlaps(hit_set.first, a_blocks, hit_set.second, false);
+ FindBlockedOverlaps(hit_set.first, a_blocks,
+ hit_set.second, false);
}
}
}
@@ -351,7 +356,8 @@ void BedIntersect::IntersectBam(string bamFile) {
hits, _sameStrand, _diffStrand,
_overlapFraction, _reciprocal);
// find the overlaps between the block in A and B
- overlapsFound = FindBlockedOverlaps(bed, bed_blocks, hits, _bamOutput);
+ overlapsFound =
+ FindBlockedOverlaps(bed, bed_blocks, hits, _bamOutput);
}
else if ((_bamOutput == false) && (_obeySplits == false))
{
diff --git a/src/utils/chromsweep/chromsweep.cpp b/src/utils/chromsweep/chromsweep.cpp
index 919e188f..88456f36 100644
--- a/src/utils/chromsweep/chromsweep.cpp
+++ b/src/utils/chromsweep/chromsweep.cpp
@@ -72,6 +72,7 @@ void ChromSweep::ScanCache() {
vector<BED>::iterator c = _cache.begin();
while (c != _cache.end())
{
+
if ((_curr_qy.chrom == c->chrom) && !(after(_curr_qy, *c))) {
if (IsValidHit(_curr_qy, *c)) {
_hits.push_back(*c);
@@ -91,9 +92,11 @@ bool ChromSweep::ChromChange()
if (_curr_qy.chrom == _curr_db.chrom) {
return false;
}
- // the query is ahead of the database. fast-forward the database to catch-up.
- else if (_curr_qy.chrom > _curr_db.chrom) {
- while (_db->GetNextBed(_curr_db, true) && _curr_db.chrom < _curr_qy.chrom)
+ // the query is ahead of the database.
+ // fast-forward the database to catch-up.
+ else if ((_curr_qy.chrom > _curr_db.chrom) && (!_db->Empty())) {
+ while (_db->GetNextBed(_curr_db, true) &&
+ _curr_db.chrom < _curr_qy.chrom)
{
}
_cache.clear();
@@ -158,7 +161,8 @@ bool ChromSweep::Next(pair<BED, vector<BED> > &next) {
if (ChromChange() == false) {
// scan the database cache for hits
ScanCache();
- // advance the db until we are ahead of the query. update hits and cache as necessary
+ // advance the db until we are ahead of the query.
+ // update hits and cache as necessary
while (!_db->Empty() &&
_curr_qy.chrom == _curr_db.chrom &&
!(after(_curr_db, _curr_qy)))
--
GitLab