Skip to content
Snippets Groups Projects
Commit 0ea1769d authored by Neil Kindlon's avatar Neil Kindlon
Browse files

Fixed bug133, missing intersection of zero length records

parent 5ba347e1
No related branches found
No related tags found
No related merge requests found
...@@ -32,7 +32,7 @@ ContextBase::ContextBase() ...@@ -32,7 +32,7 @@ ContextBase::ContextBase()
_writeOverlap(false), _writeOverlap(false),
_writeAllOverlap(false), _writeAllOverlap(false),
_haveFraction(false), _haveFraction(false),
_overlapFraction(1E-9), _overlapFraction(0.0),
_reciprocal(false), _reciprocal(false),
_sameStrand(false), _sameStrand(false),
_diffStrand(false), _diffStrand(false),
......
...@@ -160,15 +160,27 @@ bool Record::sameChromIntersects(const Record *record, ...@@ -160,15 +160,27 @@ bool Record::sameChromIntersects(const Record *record,
int otherStart = record->getStartPos(); int otherStart = record->getStartPos();
int otherEnd = record->getEndPos(); int otherEnd = record->getEndPos();
bool otherZeroLen = (otherStart - otherEnd == 0);
int maxStart = max(_startPos, otherStart); int maxStart = max(_startPos, otherStart);
int minEnd = min(_endPos, otherEnd); int minEnd = min(_endPos, otherEnd);
bool localZeroLen = (_endPos - _startPos == 0);
//rule out all cases of no intersection at all //rule out all cases of no intersection at all
if (minEnd < maxStart) { if (minEnd < maxStart) {
return false; return false;
} }
if (overlapFraction == 0.0) {
//don't care about amount of overlap.
//however, if minEnd and maxStart are equal, and
//neither record is zeroLen, return false.
if (minEnd == maxStart && !otherZeroLen && !localZeroLen) {
return false;
}
return true;
}
int overlapBases = minEnd - maxStart; int overlapBases = minEnd - maxStart;
int len = _endPos - _startPos; int len = _endPos - _startPos;
int otherLen = otherEnd - otherStart; int otherLen = otherEnd - otherStart;
......
...@@ -210,6 +210,7 @@ void RecordOutputMgr::printRecord(RecordKeyVector &keyList, RecordKeyVector *blo ...@@ -210,6 +210,7 @@ void RecordOutputMgr::printRecord(RecordKeyVector &keyList, RecordKeyVector *blo
if ((static_cast<ContextIntersect *>(_context))->getWriteAllOverlap()) { if ((static_cast<ContextIntersect *>(_context))->getWriteAllOverlap()) {
// -wao the user wants to force the reporting of 0 overlap // -wao the user wants to force the reporting of 0 overlap
if (printKeyAndTerminate(keyList)) { if (printKeyAndTerminate(keyList)) {
_currBamBlockList = NULL; _currBamBlockList = NULL;
return; return;
} }
...@@ -223,6 +224,7 @@ void RecordOutputMgr::printRecord(RecordKeyVector &keyList, RecordKeyVector *blo ...@@ -223,6 +224,7 @@ void RecordOutputMgr::printRecord(RecordKeyVector &keyList, RecordKeyVector *blo
else if ((static_cast<ContextIntersect *>(_context))->getLeftJoin()) { else if ((static_cast<ContextIntersect *>(_context))->getLeftJoin()) {
if (printKeyAndTerminate(keyList)) { if (printKeyAndTerminate(keyList)) {
_currBamBlockList = NULL; _currBamBlockList = NULL;
return; return;
} }
tab(); tab();
...@@ -230,11 +232,13 @@ void RecordOutputMgr::printRecord(RecordKeyVector &keyList, RecordKeyVector *blo ...@@ -230,11 +232,13 @@ void RecordOutputMgr::printRecord(RecordKeyVector &keyList, RecordKeyVector *blo
newline(); newline();
if (needsFlush()) flush(); if (needsFlush()) flush();
_currBamBlockList = NULL; _currBamBlockList = NULL;
return; return;
} }
} else { } else {
if (printBamRecord(keyList, true) == BAM_AS_BAM) { if (printBamRecord(keyList, true) == BAM_AS_BAM) {
_currBamBlockList = NULL; _currBamBlockList = NULL;
return; return;
} }
int hitIdx = 0; int hitIdx = 0;
...@@ -252,14 +256,17 @@ void RecordOutputMgr::printRecord(RecordKeyVector &keyList, RecordKeyVector *blo ...@@ -252,14 +256,17 @@ void RecordOutputMgr::printRecord(RecordKeyVector &keyList, RecordKeyVector *blo
newline(); newline();
} }
_currBamBlockList = NULL; _currBamBlockList = NULL;
return; return;
} else if (_context->getProgram() == ContextBase::MAP) { } else if (_context->getProgram() == ContextBase::MAP) {
printKeyAndTerminate(keyList); printKeyAndTerminate(keyList);
_currBamBlockList = NULL; _currBamBlockList = NULL;
return; return;
} else if (_context->getProgram() == ContextBase::MERGE) { } else if (_context->getProgram() == ContextBase::MERGE) {
printKeyAndTerminate(keyList); printKeyAndTerminate(keyList);
_currBamBlockList = NULL; _currBamBlockList = NULL;
return; return;
} }
} }
...@@ -380,6 +387,7 @@ void RecordOutputMgr::reportOverlapDetail(const Record *keyRecord, const Record ...@@ -380,6 +387,7 @@ void RecordOutputMgr::reportOverlapDetail(const Record *keyRecord, const Record
newline(); newline();
if (needsFlush()) flush(); if (needsFlush()) flush();
} }
const_cast<Record *>(keyRecord)->adjustZeroLength();
} }
void RecordOutputMgr::reportOverlapSummary(RecordKeyVector &keyList) void RecordOutputMgr::reportOverlapSummary(RecordKeyVector &keyList)
......
chr1 5 15 r1
chr1 7 12 r3
chr1 20 25 r2
chr1 9 9 m3
chr1 50 150 m1
chr2 20 25 m2
...@@ -764,3 +764,24 @@ $BT intersect -a oneRecordNoNewline.bed -b oneRecordNoNewline.bed > obs ...@@ -764,3 +764,24 @@ $BT intersect -a oneRecordNoNewline.bed -b oneRecordNoNewline.bed > obs
check obs exp check obs exp
rm obs exp rm obs exp
###########################################################
# Test zero length intersections in non-bam files.
############################################################
echo " intersect.new.t65...\c"
echo \
"chr1 5 15 r1 chr1 9 9 m3 0
chr1 7 12 r3 chr1 9 9 m3 0" > exp
$BT intersect -a a_testZeroLen.bed -b b_testZeroLen.bed -wo > obs
check exp obs
rm exp obs
###########################################################
# Test zero length intersections in non-bam files, -sorted
############################################################
echo " intersect.new.t66...\c"
echo \
"chr1 5 15 r1 chr1 9 9 m3 0
chr1 7 12 r3 chr1 9 9 m3 0" > exp
$BT intersect -a a_testZeroLen.bed -b b_testZeroLen.bed -wo -sorted> obs
check exp obs
rm exp obs
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment