Skip to content
Snippets Groups Projects
Commit be141f5d authored by Aaron's avatar Aaron
Browse files

Improved clarity for zero length features in coverageBed.

parent 14e45c5c
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/coverageBed/coverageBed.cpp
+ 35
− 11
View file @ be141f5d
......@@ -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);
}
}
}
......
src/utils/bedFile/bedFile.cpp
+ 30
− 6
View file @ be141f5d
......@@ -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;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment