diff --git a/src/coverageBed/coverageBed.cpp b/src/coverageBed/coverageBed.cpp
index df1a224591ced7654d0f499a2475342664d386b7..2e874845308587f47f58ca92f68c1072aa78fcd7 100644
--- a/src/coverageBed/coverageBed.cpp
+++ b/src/coverageBed/coverageBed.cpp
@@ -160,7 +160,7 @@ void BedCoverage::ReportCoverage() {
// B s ------------
int start = min(bedItr->minOverlapStart, bedItr->start);
- // track the numnber of bases in the feature covered by
+ // track the number of bases in the feature covered by
// 0, 1, 2, ... n features in A
map<unsigned int, unsigned int> depthHist;
map<unsigned int, DEPTH>::const_iterator depthItr;
@@ -181,7 +181,8 @@ void BedCoverage::ReportCoverage() {
depthHist[depth]++;
allDepthHist[depth]++;
}
- else {
+ else if ((_eachBase == true) && (bedItr->zeroLength == false))
+ {
_bedB->reportBedTab(*bedItr);
printf("%d\t%d\n", pos-bedItr->start, depth);
}
@@ -191,28 +192,51 @@ void BedCoverage::ReportCoverage() {
depth = depth - depthItr->second.ends;
}
+ // handle the special case where the user wants "per-base" depth
+ // but the current feature is length = 0.
+ if ((_eachBase == true) && (bedItr->zeroLength == true)) {
+ _bedB->reportBedTab(*bedItr);
+ printf("1\t%d\n",depth);
+ }
// Summarize the coverage for the current interval,
// assuming the user has not requested "per-base" coverage.
- if (_eachBase == false) {
+ else if (_eachBase == false)
+ {
CHRPOS length = bedItr->end - bedItr->start;
+ if (bedItr->zeroLength == true) {
+ length = 0;
+ }
totalLength += length;
int nonZeroBases = (length - zeroDepthCount);
- float fractCovered = (float) nonZeroBases / length;
-
+
+ float fractCovered = 0.0;
+ if (bedItr->zeroLength == false) {
+ fractCovered = (float) nonZeroBases / length;
+ }
+
// print a summary of the coverage
if (_writeHistogram == false) {
_bedB->reportBedTab(*bedItr);
printf("%d\t%d\t%d\t%0.7f\n", bedItr->count, nonZeroBases, length, fractCovered);
}
+ // HISTOGRAM
// report the number of bases with coverage == x
else {
- map<unsigned int, unsigned int>::const_iterator histItr = depthHist.begin();
- map<unsigned int, unsigned int>::const_iterator histEnd = depthHist.end();
- for (; histItr != histEnd; ++histItr)
- {
- float fractAtThisDepth = (float) histItr->second / length;
+ // produce a histogram when not a zero length feature.
+ if (bedItr->zeroLength == false) {
+ map<unsigned int, unsigned int>::const_iterator histItr = depthHist.begin();
+ map<unsigned int, unsigned int>::const_iterator histEnd = depthHist.end();
+ for (; histItr != histEnd; ++histItr)
+ {
+ float fractAtThisDepth = (float) histItr->second / length;
+ _bedB->reportBedTab(*bedItr);
+ printf("%d\t%d\t%d\t%0.7f\n", histItr->first, histItr->second, length, fractAtThisDepth);
+ }
+ }
+ // special case when it is a zero length feauture.
+ else {
_bedB->reportBedTab(*bedItr);
- printf("%d\t%d\t%d\t%0.7f\n", histItr->first, histItr->second, length, fractAtThisDepth);
+ printf("%d\t%d\t%d\t%0.7f\n", bedItr->count, 0, 0, 1.0000000);
}
}
}
diff --git a/src/utils/bedFile/bedFile.cpp b/src/utils/bedFile/bedFile.cpp
index 991d9015b285eeabd552610c8e25a7aca89d2f4d..914a431ade79916fcf821ed0519188b11f83703a 100644
--- a/src/utils/bedFile/bedFile.cpp
+++ b/src/utils/bedFile/bedFile.cpp
@@ -354,8 +354,15 @@ void BedFile::countHits(const BED &a, bool forceStrand) {
else if (overlaps(bedItr->start, bedItr->end, a.start, a.end) > 0) {
bedItr->count++;
- bedItr->depthMap[a.start+1].starts++;
- bedItr->depthMap[a.end].ends++;
+ if (a.zeroLength == false) {
+ bedItr->depthMap[a.start+1].starts++;
+ bedItr->depthMap[a.end].ends++;
+ }
+ else {
+ // correct for the fact that we artificially expanded the zeroLength feature
+ bedItr->depthMap[a.start+2].starts++;
+ bedItr->depthMap[a.end-1].ends++;
+ }
if (a.start < bedItr->minOverlapStart) {
bedItr->minOverlapStart = a.start;
@@ -403,8 +410,16 @@ void BedFile::countSplitHits(const vector<BED> &bedBlocks, bool forceStrand) {
continue;
}
else if (overlaps(bedItr->start, bedItr->end, blockItr->start, blockItr->end) > 0) {
- bedItr->depthMap[blockItr->start+1].starts++;
- bedItr->depthMap[blockItr->end].ends++;
+ if (blockItr->zeroLength == false) {
+ bedItr->depthMap[blockItr->start+1].starts++;
+ bedItr->depthMap[blockItr->end].ends++;
+ }
+ else {
+ // correct for the fact that we artificially expanded the zeroLength feature
+ bedItr->depthMap[blockItr->start+2].starts++;
+ bedItr->depthMap[blockItr->end-1].ends++;
+ }
+
validHits.insert(bedItr);
if (blockItr->start < bedItr->minOverlapStart)
bedItr->minOverlapStart = blockItr->start;
@@ -452,8 +467,15 @@ void BedFile::countListHits(const BED &a, int index, bool forceStrand) {
}
else if (overlaps(bedItr->start, bedItr->end, a.start, a.end) > 0) {
bedItr->counts[index]++;
- bedItr->depthMapList[index][a.start+1].starts++;
- bedItr->depthMapList[index][a.end].ends++;
+ if (a.zeroLength == false) {
+ bedItr->depthMapList[index][a.start+1].starts++;
+ bedItr->depthMapList[index][a.end].ends++;
+ }
+ else {
+ // correct for the fact that we artificially expanded the zeroLength feature
+ bedItr->depthMapList[index][a.start+2].starts++;
+ bedItr->depthMapList[index][a.end-1].ends++;
+ }
if (a.start < bedItr->minOverlapStarts[index]) {
bedItr->minOverlapStarts[index] = a.start;
@@ -522,6 +544,7 @@ void BedFile::loadBedCovFileIntoMap() {
bedCov.score = bedEntry.score;
bedCov.strand = bedEntry.strand;
bedCov.otherFields = bedEntry.otherFields;
+ bedCov.zeroLength = bedEntry.zeroLength;
bedCov.count = 0;
bedCov.minOverlapStart = INT_MAX;
@@ -551,6 +574,7 @@ void BedFile::loadBedCovListFileIntoMap() {
bedCovList.score = bedEntry.score;
bedCovList.strand = bedEntry.strand;
bedCovList.otherFields = bedEntry.otherFields;
+ bedCovList.zeroLength = bedEntry.zeroLength;
bedCovListMap[bedEntry.chrom][bin].push_back(bedCovList);
bedEntry = nullBed;