diff --git a/src/bamToBed/bamToBed.cpp b/src/bamToBed/bamToBed.cpp
index 344be8622634eec92fad740a74b4664dc6880a22..bcf7c2b4246a7d807c20645858ad9794d45ce30f 100644
--- a/src/bamToBed/bamToBed.cpp
+++ b/src/bamToBed/bamToBed.cpp
@@ -413,7 +413,7 @@ void PrintBed(const BamAlignment &bam, const RefVector &refs, bool useEditDista
vector<BED> bedBlocks;
string chrom = refs.at(bam.RefID).RefName;
// extract the block starts and lengths from the CIGAR string
- GetBamBlocks(bam, chrom, bedBlocks);
+ GetBamBlocks(bam, chrom, bedBlocks, false, true);
unsigned int i;
for (i = 0; i < bedBlocks.size(); ++i) {
@@ -442,7 +442,7 @@ void PrintBed12(const BamAlignment &bam, const RefVector &refs, bool useEditDist
string chrom = refs.at(bam.RefID).RefName;
CHRPOS alignmentEnd = bam.GetEndPosition();
// extract the block starts and lengths from the CIGAR string
- GetBamBlocks(bam, chrom, bedBlocks);
+ GetBamBlocks(bam, chrom, bedBlocks, false, true);
// write BED6 portion
if (useEditDistance == false && bamTag == "") {
diff --git a/src/coverageBed/coverageBed.cpp b/src/coverageBed/coverageBed.cpp
index dfd3a199ddc84be80a4a22de6ce1428f2e95827e..2d3bba81c873aa28a99609dc4eac3ecdd40cf6ad 100644
--- a/src/coverageBed/coverageBed.cpp
+++ b/src/coverageBed/coverageBed.cpp
@@ -117,7 +117,7 @@ void BedCoverage::CollectCoverageBam(string bamFile) {
bedVector bedBlocks;
// since we are counting coverage, we do want to split blocks when a
// deletion (D) CIGAR op is encountered (hence the true for the last parm)
- GetBamBlocks(bam, refs.at(bam.RefID).RefName, bedBlocks, true);
+ GetBamBlocks(bam, refs.at(bam.RefID).RefName, bedBlocks, false, true);
// use countSplitHits to avoid over-counting each split chunk
// as distinct read coverage.
_bedB->countSplitHits(bedBlocks, _sameStrand, _diffStrand, _countsOnly);
diff --git a/src/genomeCoverageBed/genomeCoverageBed.cpp b/src/genomeCoverageBed/genomeCoverageBed.cpp
index a06ac20425f07b5fe4d254cdebcdf1c96229ae17..cdb044464ec49e96c4cd2b4b4d256a41e8354baa 100644
--- a/src/genomeCoverageBed/genomeCoverageBed.cpp
+++ b/src/genomeCoverageBed/genomeCoverageBed.cpp
@@ -243,11 +243,16 @@ void BedGenomeCoverage::CoverageBam(string bamFile) {
StartNewChrom(chrom);
// add coverage accordingly.
- if (_obeySplits) {
+ if (!_only_5p_end && !_only_3p_end) {
bedVector bedBlocks;
- // since we are counting coverage, we do want to split blocks when a
- // deletion (D) CIGAR op is encountered (hence the true for the last parm)
- GetBamBlocks(bam, refs.at(bam.RefID).RefName, bedBlocks, true);
+ // we always want to split blocks when a D CIGAR op is found.
+ // if the user invokes -split, we want to also split on N ops.
+ if (_obeySplits) { // "D" true, "N" true
+ GetBamBlocks(bam, refs.at(bam.RefID).RefName, bedBlocks, true, true);
+ }
+ else { // "D" true, "N" false
+ GetBamBlocks(bam, refs.at(bam.RefID).RefName, bedBlocks, true, false);
+ }
AddBlockedCoverage(bedBlocks);
}
else if (_only_5p_end) {
@@ -258,8 +263,6 @@ void BedGenomeCoverage::CoverageBam(string bamFile) {
int pos = ( bam.IsReverseStrand() ) ? start : end;
AddCoverage(pos,pos);
}
- else
- AddCoverage(start, end);
}
// close the BAM
reader.Close();
diff --git a/src/intersectBed/intersectBed.cpp b/src/intersectBed/intersectBed.cpp
index 559eec23c4fb5659fd08b5b26bb98064d84fbcc9..1d410ae16f23994ef30ee00ac68b24fa4f232bcb 100644
--- a/src/intersectBed/intersectBed.cpp
+++ b/src/intersectBed/intersectBed.cpp
@@ -323,7 +323,7 @@ void BedIntersect::IntersectBam(string bamFile) {
// break alignment into discrete blocks,
bedVector bed_blocks;
string chrom = refs.at(bam.RefID).RefName;
- GetBamBlocks(bam, chrom, bed_blocks);
+ GetBamBlocks(bam, chrom, bed_blocks, false, true);
// create a basic BED entry from the BAM alignment
BED bed;
MakeBedFromBam(bam, chrom, bed_blocks, bed);
diff --git a/src/utils/BlockedIntervals/BlockedIntervals.cpp b/src/utils/BlockedIntervals/BlockedIntervals.cpp
index eead4ca77659e2d52d0586c2ad6044120fa94c6a..659f299b046291e282bec7f9c6a9722c35e56fb5 100644
--- a/src/utils/BlockedIntervals/BlockedIntervals.cpp
+++ b/src/utils/BlockedIntervals/BlockedIntervals.cpp
@@ -15,7 +15,8 @@
void GetBamBlocks(const BamAlignment &bam,
const string &chrom,
bedVector &bedBlocks,
- bool breakOnDeletionOps)
+ bool breakOnDeletionOps,
+ bool breakOnSkipOps)
{
vector<int> starts;
vector<int> lengths;
@@ -47,10 +48,14 @@ void GetBamBlocks(const BamAlignment &bam,
}
case ('P') : break;
case ('N') :
- bedBlocks.push_back( BED(chrom, currPosition, currPosition + blockLength,
- bam.Name, ToString(bam.MapQuality), strand) );
- currPosition += cigItr->Length + blockLength;
- blockLength = 0;
+ if (!breakOnSkipOps)
+ blockLength += cigItr->Length;
+ else {
+ bedBlocks.push_back( BED(chrom, currPosition, currPosition + blockLength,
+ bam.Name, ToString(bam.MapQuality), strand) );
+ currPosition += cigItr->Length + blockLength;
+ blockLength = 0;
+ }
case ('H') : break; // for 'H' - do nothing, move to next op
default :
printf("ERROR: Invalid Cigar op type\n"); // shouldn't get here
diff --git a/src/utils/BlockedIntervals/BlockedIntervals.h b/src/utils/BlockedIntervals/BlockedIntervals.h
index 5a3a5c2cbfea4c5a39e2874171c564650904a2ab..1f6c182e3686d8f44bb87f2d4f1caed9aef5d1c4 100644
--- a/src/utils/BlockedIntervals/BlockedIntervals.h
+++ b/src/utils/BlockedIntervals/BlockedIntervals.h
@@ -21,9 +21,10 @@ using namespace BamTools;
into discrete alignment blocks.
*/
void GetBamBlocks(const BamAlignment &bam,
- const string &chrom,
- bedVector &bedBlocks,
- bool breakOnDeletionOps = false);
+ const string &chrom,
+ bedVector &bedBlocks,
+ bool breakOnDeletionOps,
+ bool breakOnSkipOps);
/* break a BED12 record into discrete BED6 blocks. */
void GetBedBlocks(const BED &bed, bedVector &bedBlocks);