Commit 2f6dfe9c authored by Aaron's avatar Aaron
Browse files

First cut at BlockedIntervals class to cleanup spliced BAMs and BED12

parent b4a33ce9
......@@ -61,6 +61,7 @@ UTIL_SUBDIRS = $(SRC_DIR)/utils/lineFileUtilities \
$(SRC_DIR)/utils/tabFile \
$(SRC_DIR)/utils/BamTools \
$(SRC_DIR)/utils/BamTools-Ancillary \
$(SRC_DIR)/utils/BlockedIntervals \
$(SRC_DIR)/utils/Fasta \
$(SRC_DIR)/utils/VectorOps \
$(SRC_DIR)/utils/genomeFile
......
......@@ -10,7 +10,7 @@ INCLUDES = -I$(UTILITIES_DIR)/bedFile/ \
-I$(UTILITIES_DIR)/lineFileUtilities/ \
-I$(UTILITIES_DIR)/fileType/ \
-I$(UTILITIES_DIR)/BamTools/include \
-I$(UTILITIES_DIR)/BamTools-Ancillary \
-I$(UTILITIES_DIR)/BlockedIntervals \
-I$(UTILITIES_DIR)/version/
......@@ -19,7 +19,7 @@ INCLUDES = -I$(UTILITIES_DIR)/bedFile/ \
# ----------------------------------
SOURCES= bamToBed.cpp
OBJECTS= $(SOURCES:.cpp=.o)
_EXT_OBJECTS=BamAncillary.o
_EXT_OBJECTS=BamAncillary.o BlockedIntervals.o
EXT_OBJECTS=$(patsubst %,$(OBJ_DIR)/%,$(_EXT_OBJECTS))
BUILT_OBJECTS= $(patsubst %,$(OBJ_DIR)/%,$(OBJECTS))
PROGRAM= bamToBed
......@@ -35,6 +35,7 @@ $(BUILT_OBJECTS): $(SOURCES)
$(EXT_OBJECTS):
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/BamAncillary/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/BlockedIntervals/
clean:
@echo "Cleaning up."
......
......@@ -11,7 +11,7 @@
******************************************************************************/
#include "api/BamReader.h"
#include "api/BamAux.h"
#include "BamAncillary.h"
#include "BlockedIntervals.h"
#include "bedFile.h"
#include "version.h"
using namespace BamTools;
......@@ -289,50 +289,10 @@ void ConvertBamToBedpe(const string &bamFile, const bool &useEditDistance) {
}
void ParseCigarBed12(const vector<CigarOp> &cigar, vector<int> &blockStarts, vector<int> &blockLengths, unsigned int &alignmentEnd) {
int currPosition = 0;
int blockLength = 0;
// Rip through the CIGAR ops and figure out if there is more
// than one block for this alignment
vector<CigarOp>::const_iterator cigItr = cigar.begin();
vector<CigarOp>::const_iterator cigEnd = cigar.end();
for (; cigItr != cigEnd; ++cigItr) {
switch (cigItr->Type) {
case ('M') :
blockLength += cigItr->Length;
currPosition += cigItr->Length;
case ('I') : break;
case ('S') : break;
case ('D') : break;
blockLength += cigItr->Length;
currPosition += cigItr->Length;
case ('P') : break;
case ('N') :
blockStarts.push_back(currPosition + cigItr->Length);
blockLengths.push_back(blockLength);
currPosition += cigItr->Length;
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
exit(1);
}
}
// add the kast block and set the
// alignment end (i.e., relative to the start)
blockLengths.push_back(blockLength);
alignmentEnd = currPosition;
}
string BuildCigarString(const vector<CigarOp> &cigar) {
stringstream cigarString;
for (size_t i = 0; i < cigar.size(); ++i) {
//cerr << cigar[i].Type << " " << cigar[i].Length << endl;
switch (cigar[i].Type) {
case ('M') :
case ('I') :
......@@ -449,16 +409,18 @@ void PrintBed(const BamAlignment &bam, const RefVector &refs, bool useEditDista
// Report each chunk of the BAM alignment as a discrete BED entry
// For example 10M100N10M would be reported as two seprate BED entries of length 10
else {
// parse the CIGAR string and figure out the alignment blocks
vector<BED> bedBlocks;
// Load the alignment blocks in bam into the bedBlocks vector.
// Don't trigger a new block when a "D" (deletion) CIGAR op is found.
getBamBlocks(bam, refs, bedBlocks, false);
vector<BED>::const_iterator bedItr = bedBlocks.begin();
vector<BED>::const_iterator bedEnd = bedBlocks.end();
for (; bedItr != bedEnd; ++bedItr) {
printf("%s\t%d\t%d\t\%s\t%d\t%s\n", refs.at(bam.RefID).RefName.c_str(), bedItr->start,
bedItr->end, name.c_str(), bam.MapQuality, strand.c_str());
string chrom = refs.at(bam.RefID).RefName;
// extract the block starts and lengths from the CIGAR string
GetBamBlocks(bam, chrom, bedBlocks);
unsigned int i;
for (i = 0; i < bedBlocks.size(); ++i) {
BED curr = bedBlocks[i];
printf("%s\t%d\t%d\t\%s\t%d\t%s\n",
chrom.c_str(), curr.start, curr.end,
name.c_str(), bam.MapQuality, strand.c_str());
}
}
}
......@@ -476,14 +438,11 @@ void PrintBed12(const BamAlignment &bam, const RefVector &refs, bool useEditDist
if (bam.IsSecondMate()) name += "/2";
// parse the CIGAR string and figure out the alignment blocks
unsigned int alignmentEnd;
vector<int> blockLengths;
vector<int> blockStarts;
blockStarts.push_back(0);
vector<BED> bedBlocks;
string chrom = refs.at(bam.RefID).RefName;
CHRPOS alignmentEnd = bam.GetEndPosition();
// extract the block starts and lengths from the CIGAR string
ParseCigarBed12(bam.CigarData, blockStarts, blockLengths, alignmentEnd);
alignmentEnd += bam.Position;
GetBamBlocks(bam, chrom, bedBlocks);
// write BED6 portion
if (useEditDistance == false && bamTag == "") {
......@@ -514,20 +473,20 @@ void PrintBed12(const BamAlignment &bam, const RefVector &refs, bool useEditDist
}
// write the colors, etc.
printf("%d\t%d\t%s\t%d\t", bam.Position, alignmentEnd, color.c_str(), (int) blockStarts.size());
printf("%d\t%d\t%s\t%d\t", bam.Position, alignmentEnd, color.c_str(), (int) bedBlocks.size());
// now write the lengths portion
unsigned int b;
for (b = 0; b < blockLengths.size() - 1; ++b) {
printf("%d,", blockLengths[b]);
for (b = 0; b < bedBlocks.size() - 1; ++b) {
printf("%d,", bedBlocks[b].end - bedBlocks[b].start);
}
printf("%d\t", blockLengths[b]);
printf("%d\t", bedBlocks[b].end - bedBlocks[b].start);
// now write the starts portion
for (b = 0; b < blockStarts.size() - 1; ++b) {
printf("%d,", blockStarts[b]);
for (b = 0; b < bedBlocks.size() - 1; ++b) {
printf("%d,", bedBlocks[b].start - bam.Position);
}
printf("%d\n", blockStarts[b]);
printf("%d\n", bedBlocks[b].start - bam.Position);
}
......
......@@ -9,6 +9,8 @@ INCLUDES = -I$(UTILITIES_DIR)/bedFile/ \
-I$(UTILITIES_DIR)/lineFileUtilities/ \
-I$(UTILITIES_DIR)/gzstream/ \
-I$(UTILITIES_DIR)/fileType/ \
-I$(UTILITIES_DIR)/BamTools/include \
-I$(UTILITIES_DIR)/BlockedIntervals \
-I$(UTILITIES_DIR)/version/
# ----------------------------------
......@@ -16,7 +18,7 @@ INCLUDES = -I$(UTILITIES_DIR)/bedFile/ \
# ----------------------------------
SOURCES= bed12ToBed6.cpp
OBJECTS= $(SOURCES:.cpp=.o)
_EXT_OBJECTS=bedFile.o lineFileUtilities.o gzstream.o fileType.o
_EXT_OBJECTS=bedFile.o BlockedIntervals.o lineFileUtilities.o gzstream.o fileType.o
EXT_OBJECTS=$(patsubst %,$(OBJ_DIR)/%,$(_EXT_OBJECTS))
BUILT_OBJECTS= $(patsubst %,$(OBJ_DIR)/%,$(OBJECTS))
PROGRAM= bed12ToBed6
......@@ -35,6 +37,7 @@ $(EXT_OBJECTS):
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/bedFile/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/gzstream/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/fileType/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/BlockedIntervals/
clean:
@echo "Cleaning up."
......
......@@ -10,6 +10,7 @@
Licenced under the GNU General Public License 2.0 license.
******************************************************************************/
#include "lineFileUtilities.h"
#include "BlockedIntervals.h"
#include "bedFile.h"
#include "version.h"
#include <vector>
......@@ -140,7 +141,7 @@ void ProcessBed(istream &bedInput, BedFile *bed) {
if (bed->_status == BED_VALID) {
bedVector bedBlocks; // vec to store the discrete BED "blocks" from a
splitBedIntoBlocks(bedEntry, bedBlocks);
GetBedBlocks(bedEntry, bedBlocks);
for (int i = 0; i < (int) bedBlocks.size(); ++i) {
if (addBlockNums == false) {
......
......@@ -11,15 +11,15 @@ INCLUDES = -I$(UTILITIES_DIR)/bedFile/ \
-I$(UTILITIES_DIR)/lineFileUtilities/ \
-I$(UTILITIES_DIR)/fileType/ \
-I$(UTILITIES_DIR)/BamTools/include \
-I$(UTILITIES_DIR)/BamTools-Ancillary \
-I$(UTILITIES_DIR)/BlockedIntervals \
-I$(UTILITIES_DIR)/version/
# ----------------------------------
# define our source and object files
# ----------------------------------
SOURCES= coverageMain.cpp coverageBed.cpp
OBJECTS= $(SOURCES:.cpp=.o)
_EXT_OBJECTS=bedFile.o lineFileUtilities.o gzstream.o fileType.o BamAncillary.o
EXT_OBJECTS=$(patsubst %,$(OBJ_DIR)/%,$(_EXT_OBJECTS))
# _EXT_OBJECTS=bedFile.o lineFileUtilities.o gzstream.o fileType.o BlockedIntervals.o
# EXT_OBJECTS=$(patsubst %,$(OBJ_DIR)/%,$(_EXT_OBJECTS))
BUILT_OBJECTS= $(patsubst %,$(OBJ_DIR)/%,$(OBJECTS))
PROGRAM= coverageBed
......@@ -32,13 +32,13 @@ $(BUILT_OBJECTS): $(SOURCES)
@echo " * compiling" $(*F).cpp
@$(CXX) -c -o $@ $(*F).cpp $(LDFLAGS) $(CXXFLAGS) $(INCLUDES)
$(EXT_OBJECTS):
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/bedFile/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/lineFileUtilities/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/BamTools/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/BamTools-Ancillary/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/gzstream/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/fileType/
# $(EXT_OBJECTS):
# @$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/bedFile/
# @$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/lineFileUtilities/
# @$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/BamTools/
# @$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/BlockedIntervals/
# @$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/gzstream/
# @$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/fileType/
clean:
@echo "Cleaning up."
......
......@@ -62,7 +62,7 @@ void BedCoverage::CollectCoverageBed() {
// split the BED into discrete blocksand process each independently.
else {
bedVector bedBlocks;
splitBedIntoBlocks(a, bedBlocks);
GetBedBlocks(a, bedBlocks);
// use countSplitHits to avoid over-counting each split chunk
// as distinct read coverage.
_bedB->countSplitHits(bedBlocks, _sameStrand, _diffStrand, _countsOnly);
......@@ -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, bedBlocks, true);
GetBamBlocks(bam, refs.at(bam.RefID).RefName, bedBlocks, true);
// use countSplitHits to avoid over-counting each split chunk
// as distinct read coverage.
_bedB->countSplitHits(bedBlocks, _sameStrand, _diffStrand, _countsOnly);
......
......@@ -16,7 +16,7 @@
#include "api/BamReader.h"
#include "api/BamAux.h"
#include "BamAncillary.h"
#include "BlockedIntervals.h"
using namespace BamTools;
#include <vector>
......
......@@ -11,7 +11,7 @@ INCLUDES = -I$(UTILITIES_DIR)/bedFile/ \
-I$(UTILITIES_DIR)/lineFileUtilities/ \
-I$(UTILITIES_DIR)/fileType/ \
-I$(UTILITIES_DIR)/BamTools/include \
-I$(UTILITIES_DIR)/BamTools-Ancillary \
-I$(UTILITIES_DIR)/BlockedIntervals \
-I$(UTILITIES_DIR)/version/
# ----------------------------------
# define our source and object files
......@@ -37,7 +37,7 @@ $(EXT_OBJECTS):
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/bedFile/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/genomeFile
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/BamTools/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/BamTools-Ancillary/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/BlockedIntervals/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/gzstream/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/fileType/
......
......@@ -172,7 +172,7 @@ void BedGenomeCoverage::CoverageBed() {
if (_obeySplits == true) {
bedVector bedBlocks; // vec to store the discrete BED "blocks"
splitBedIntoBlocks(a, bedBlocks);
GetBedBlocks(a, bedBlocks);
AddBlockedCoverage(bedBlocks);
}
else if (_only_5p_end) {
......@@ -247,7 +247,7 @@ void BedGenomeCoverage::CoverageBam(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, bedBlocks, true);
GetBamBlocks(bam, refs.at(bam.RefID).RefName, bedBlocks, true);
AddBlockedCoverage(bedBlocks);
}
else if (_only_5p_end) {
......
......@@ -12,7 +12,7 @@ Licenced under the GNU General Public License 2.0 license.
#include "bedFile.h"
#include "genomeFile.h"
#include "BamAncillary.h"
#include "BlockedIntervals.h"
#include "api/BamReader.h"
#include "api/BamAux.h"
using namespace BamTools;
......
......@@ -11,7 +11,7 @@ INCLUDES = -I$(UTILITIES_DIR)/bedFile/ \
-I$(UTILITIES_DIR)/lineFileUtilities/ \
-I$(UTILITIES_DIR)/fileType/ \
-I$(UTILITIES_DIR)/BamTools/include \
-I$(UTILITIES_DIR)/BamTools-Ancillary \
-I$(UTILITIES_DIR)/BlockedIntervals \
-I$(UTILITIES_DIR)/chromsweep \
-I$(UTILITIES_DIR)/version/
......@@ -20,7 +20,7 @@ INCLUDES = -I$(UTILITIES_DIR)/bedFile/ \
# ----------------------------------
SOURCES= intersectMain.cpp intersectBed.cpp
OBJECTS= $(SOURCES:.cpp=.o)
_EXT_OBJECTS=bedFile.o lineFileUtilities.o BamAncillary.o gzstream.o fileType.o chromsweep.o
_EXT_OBJECTS=bedFile.o lineFileUtilities.o BlockedIntervals.o gzstream.o fileType.o chromsweep.o
EXT_OBJECTS=$(patsubst %,$(OBJ_DIR)/%,$(_EXT_OBJECTS))
BUILT_OBJECTS= $(patsubst %,$(OBJ_DIR)/%,$(OBJECTS))
PROGRAM= intersectBed
......@@ -37,7 +37,7 @@ $(EXT_OBJECTS):
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/bedFile/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/lineFileUtilities/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/BamTools/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/BamTools-Ancillary/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/BlockedIntervals/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/gzstream/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/fileType/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/chromsweep/
......
......@@ -17,47 +17,6 @@ Helper functions
************************************/
bool BedIntersect::processHits(const BED &a, const vector<BED> &hits) {
// // how many overlaps are there b/w the bed and the set of hits?
// CHRPOS s, e;
// int overlapBases;
// int numOverlaps = 0;
// bool hitsFound = false;
// int aLength = (a.end - a.start); // the length of a in b.p.
//
// // loop through the hits and report those that meet the user's criteria
// vector<BED>::const_iterator h = hits.begin();
// vector<BED>::const_iterator hitsEnd = hits.end();
// for (; h != hitsEnd; ++h) {
// s = max(a.start, h->start);
// e = min(a.end, h->end);
// overlapBases = (e - s); // the number of overlapping bases b/w a and b
//
// // is there enough overlap relative to the user's request? (default ~ 1bp)
// if ( ( (float) overlapBases / (float) aLength ) >= _overlapFraction ) {
// // Report the hit if the user doesn't care about reciprocal overlap between A and B.
// if (_reciprocal == false) {
// hitsFound = true;
// numOverlaps++;
// if (_printable == true)
// ReportOverlapDetail(overlapBases, a, *h, s, e);
// }
// // we require there to be sufficient __reciprocal__ overlap
// else {
// int bLength = (h->end - h->start);
// float bOverlap = ( (float) overlapBases / (float) bLength );
// if (bOverlap >= _overlapFraction) {
// hitsFound = true;
// numOverlaps++;
// if (_printable == true)
// ReportOverlapDetail(overlapBases, a, *h, s, e);
// }
// }
// }
// }
// // report the summary of the overlaps if requested.
// ReportOverlapSummary(a, numOverlaps);
// // were hits found for this BED feature?
// return hitsFound;
bool hitsFound = false;
if (_printable == true) {
vector<BED>::const_iterator h = hits.begin();
......@@ -220,7 +179,7 @@ void BedIntersect::IntersectBed() {
// split the BED12 into blocks and look for overlaps in each discrete block
else {
bedVector bedBlocks; // vec to store the discrete BED "blocks"
splitBedIntoBlocks(a,bedBlocks);
GetBedBlocks(a,bedBlocks);
vector<BED>::const_iterator bedItr = bedBlocks.begin();
vector<BED>::const_iterator bedEnd = bedBlocks.end();
......@@ -322,7 +281,7 @@ void BedIntersect::IntersectBam(string bamFile) {
bool overlapFoundForBlock;
bedVector bedBlocks; // vec to store the discrete BED "blocks" from a
// we don't want to split on "D" ops, hence the "false"
getBamBlocks(bam, refs, bedBlocks, false);
GetBamBlocks(bam, refs.at(bam.RefID).RefName, bedBlocks);
vector<BED>::const_iterator bedItr = bedBlocks.begin();
vector<BED>::const_iterator bedEnd = bedBlocks.end();
......@@ -354,7 +313,7 @@ void BedIntersect::IntersectBam(string bamFile) {
// look for overlaps only within each block.
else {
bedVector bedBlocks; // vec to store the discrete BED "blocks" from a
getBamBlocks(bam, refs, bedBlocks, false);
GetBamBlocks(bam, refs.at(bam.RefID).RefName, bedBlocks);
vector<BED>::const_iterator bedItr = bedBlocks.begin();
vector<BED>::const_iterator bedEnd = bedBlocks.end();
......
......@@ -17,7 +17,7 @@
#include "api/BamReader.h"
#include "api/BamWriter.h"
#include "api/BamAux.h"
#include "BamAncillary.h"
#include "BlockedIntervals.h"
using namespace BamTools;
......
/*****************************************************************************
BlockedIntervals.h
(c) 2009 - Aaron Quinlan
Hall Laboratory
Department of Biochemistry and Molecular Genetics
University of Virginia
aaronquinlan@gmail.com
Licensed under the GNU General Public License 2.0 license.
******************************************************************************/
#include "BlockedIntervals.h"
void GetBamBlocks(const BamAlignment &bam,
const string &chrom,
vector<BED> &bedBlocks,
bool breakOnDeletionOps)
{
vector<int> starts;
vector<int> lengths;
starts.push_back(0);
string strand;
bam.IsReverseStrand() ? strand = "-" : strand = "+";
CHRPOS currPosition = bam.Position;
int blockLength = 0;
// Rip through the CIGAR ops and figure out if there is more
// than one block for this alignment
vector<CigarOp>::const_iterator cigItr = bam.CigarData.begin();
vector<CigarOp>::const_iterator cigEnd = bam.CigarData.end();
for (; cigItr != cigEnd; ++cigItr) {
switch (cigItr->Type) {
case ('M') :
blockLength += cigItr->Length;
case ('I') : break;
case ('S') : break;
case ('D') :
if (!breakOnDeletionOps)
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 ('P') : break;
case ('N') :
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
exit(1);
}
}
bedBlocks.push_back( BED(chrom, currPosition, currPosition + blockLength,
bam.Name, ToString(bam.MapQuality), strand) );
}
void GetBedBlocks(const BED &bed, bedVector &bedBlocks) {
if (bed.fields.size() < 6) {
cerr << "Input error: Cannot split into blocks. Found interval with fewer than 12 columns." << endl;
exit(1);
}
int blockCount = atoi(bed.fields[9].c_str());
if ( blockCount <= 0 ) {
cerr << "Input error: found interval having <= 0 blocks." << endl;
exit(1);
}
else if ( blockCount == 1 ) {
//take a short-cut for single blocks
bedBlocks.push_back(bed);
}
else {
// get the comma-delimited strings for the BED12 block starts and block ends.
string blockSizes(bed.fields[10]);
string blockStarts(bed.fields[11]);
vector<int> sizes;
vector<int> starts;
Tokenize(blockSizes, sizes, ",");
Tokenize(blockStarts, starts, ",");
if ( sizes.size() != (size_t) blockCount || starts.size() != (size_t) blockCount ) {
cerr << "Input error: found interval with block-counts not matching starts/sizes on line." << endl;
exit(1);
}
// add each BED block to the bedBlocks vector
for (UINT i = 0; i < (UINT) blockCount; ++i) {
CHRPOS blockStart = bed.start + starts[i];
CHRPOS blockEnd = bed.start + starts[i] + sizes[i];
BED currBedBlock(bed.chrom, blockStart, blockEnd,
bed.name, bed.score, bed.strand, bed.fields, bed.other_idxs);
bedBlocks.push_back(currBedBlock);
}
}
}
\ No newline at end of file
/*****************************************************************************
BlockedIntervals.h
(c) 2009 - Aaron Quinlan
Hall Laboratory
Department of Biochemistry and Molecular Genetics
University of Virginia
aaronquinlan@gmail.com
Licensed under the GNU General Public License 2.0 license.
******************************************************************************/
#include <vector>
#include "bedFile.h"
#include "api/BamAlignment.h"
using namespace std;
using namespace BamTools;
/*
Using the CIGAR string, break a BAM alignment
into discrete alignment blocks.
*/
void GetBamBlocks(const BamAlignment &bam,
const string &chrom,
vector<BED> &bedBlocks,
bool breakOnDeletionOps = false);
/* break a BED12 record into discrete BED6 blocks. */
void GetBedBlocks(const BED &bed, bedVector &bedBlocks);
\ No newline at end of file
OBJ_DIR = ../../../obj/
BIN_DIR = ../../../bin/
UTILITIES_DIR = ../
INCLUDES = -I$(UTILITIES_DIR)/BamTools/include \
-I$(UTILITIES_DIR)/bedFile/ \
-I$(UTILITIES_DIR)/gzstream/ \
-I$(UTILITIES_DIR)/fileType/ \
-I$(UTILITIES_DIR)/lineFileUtilities/
# ----------------------------------
# define our source and object files
# ----------------------------------
SOURCES= BlockedIntervals.cpp
OBJECTS= $(SOURCES:.cpp=.o)
_EXT_OBJECTS=bedFile.o fileType.o gzstream.o
EXT_OBJECTS=$(patsubst %,$(OBJ_DIR)/%,$(_EXT_OBJECTS))
BUILT_OBJECTS= $(patsubst %,$(OBJ_DIR)/%,$(OBJECTS))
all: $(BUILT_OBJECTS)
.PHONY: all
$(BUILT_OBJECTS): $(SOURCES)
@echo " * compiling" $(*F).cpp
@$(CXX) -c -o $@ $(*F).cpp $(LDFLAGS) $(CXXFLAGS) $(INCLUDES) -L$(BT_ROOT)/lib
$(EXT_OBJECTS):
@$(MAKE) --no-print-directory -C $(INCLUDES)