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

Use new GetNextBed() interface.

parent dfb191ce
No related branches found
No related tags found
No related merge requests found
UTILITIES_DIR = ../utils/
OBJ_DIR = ../../obj/
BIN_DIR = ../../bin/
# -------------------
# define our includes
# -------------------
INCLUDES = -I$(UTILITIES_DIR)/bedFile/ -I$(UTILITIES_DIR)/sequenceUtilities/ -I$(UTILITIES_DIR)/lineFileUtilities/ -I$(UTILITIES_DIR)/version/ -I$(UTILITIES_DIR)/gzstream/ -I$(UTILITIES_DIR)/fileType/
# ----------------------------------
# define our source and object files
# ----------------------------------
SOURCES= cuffToTransMain.cpp cuffToTrans.cpp Fasta.cpp split.cpp
OBJECTS= $(SOURCES:.cpp=.o)
_EXT_OBJECTS=bedFile.o sequenceUtils.o lineFileUtilities.o gzstream.o fileType.o
EXT_OBJECTS=$(patsubst %,$(OBJ_DIR)/%,$(_EXT_OBJECTS))
BUILT_OBJECTS= $(patsubst %,$(OBJ_DIR)/%,$(OBJECTS))
PROGRAM= cuffToTrans
all: $(PROGRAM)
.PHONY: all
$(PROGRAM): $(BUILT_OBJECTS) $(EXT_OBJECTS)
@echo " * linking $(PROGRAM)"
@$(CXX) $(LDFLAGS) $(CXXFLAGS) -o $(BIN_DIR)/$@ $^ $(LIBS)
$(BUILT_OBJECTS): $(SOURCES)
@echo " * compiling" $(*F).cpp
@$(CXX) -c -o $@ $(*F).cpp $(LDFLAGS) $(CXXFLAGS) $(INCLUDES)
$(EXT_OBJECTS):
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/sequenceUtilities/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/lineFileUtilities/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/bedFile/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/gzstream/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/fileType/
clean:
@echo "Cleaning up."
@rm -f $(OBJ_DIR)/* $(BIN_DIR)/*
.PHONY: clean
\ No newline at end of file
UTILITIES_DIR = ../utils/
OBJ_DIR = ../../obj/
BIN_DIR = ../../bin/
# -------------------
# define our includes
# -------------------
INCLUDES = -I$(UTILITIES_DIR)/bedFile/ -I$(UTILITIES_DIR)/lineFileUtilities/ -I$(UTILITIES_DIR)/version/ -I$(UTILITIES_DIR)/gzstream/ -I$(UTILITIES_DIR)/fileType/
# ----------------------------------
# define our source and object files
# ----------------------------------
SOURCES= fjoinMain.cpp fjoin.cpp
OBJECTS= $(SOURCES:.cpp=.o)
_EXT_OBJECTS=bedFile.o lineFileUtilities.o gzstream.o fileType.o
EXT_OBJECTS=$(patsubst %,$(OBJ_DIR)/%,$(_EXT_OBJECTS))
BUILT_OBJECTS= $(patsubst %,$(OBJ_DIR)/%,$(OBJECTS))
PROGRAM= fjoin
all: $(PROGRAM)
.PHONY: all
$(PROGRAM): $(BUILT_OBJECTS) $(EXT_OBJECTS)
@echo " * linking $(PROGRAM)"
@$(CXX) $(LDFLAGS) $(CXXFLAGS) -o $(BIN_DIR)/$@ $^ $(LIBS)
$(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)/gzstream/
@$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/fileType/
clean:
@echo "Cleaning up."
@rm -f $(OBJ_DIR)/* $(BIN_DIR)/*
.PHONY: clean
/*****************************************************************************
intersectBed.cpp
(c) 2009 - Aaron Quinlan
Hall Laboratory
Department of Biochemistry and Molecular Genetics
University of Virginia
aaronquinlan@gmail.com
Licenced under the GNU General Public License 2.0 license.
******************************************************************************/
#include "lineFileUtilities.h"
#include "fjoin.h"
#include <queue>
#include <set>
bool leftOf(const BED &a, const BED &b);
bool BedIntersect::processHits(BED &a, vector<BED> &hits) {
// how many overlaps are there b/w the bed and the set of hits?
int s, e, 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;
}
/*
Constructor
*/
BedIntersect::BedIntersect(string bedAFile, string bedBFile, bool anyHit,
bool writeA, bool writeB, bool writeOverlap, bool writeAllOverlap,
float overlapFraction, bool noHit, bool writeCount, bool forceStrand,
bool reciprocal, bool obeySplits, bool bamInput, bool bamOutput) {
_bedAFile = bedAFile;
_bedBFile = bedBFile;
_anyHit = anyHit;
_noHit = noHit;
_writeA = writeA;
_writeB = writeB;
_writeOverlap = writeOverlap;
_writeAllOverlap = writeAllOverlap;
_writeCount = writeCount;
_overlapFraction = overlapFraction;
_forceStrand = forceStrand;
_reciprocal = reciprocal;
_obeySplits = obeySplits;
_bamInput = bamInput;
_bamOutput = bamOutput;
if (_anyHit || _noHit || _writeCount)
_printable = false;
else
_printable = true;
// create new BED file objects for A and B
_bedA = new BedFile(bedAFile);
_bedB = new BedFile(bedBFile);
IntersectBed();
}
/*
Destructor
*/
BedIntersect::~BedIntersect(void) {
}
bool leftOf(const BED &a, const BED &b) {
return (a.end <= b.start);
}
void BedIntersect::ReportOverlapDetail(const int &overlapBases, const BED &a, const BED &b,
const CHRPOS &s, const CHRPOS &e) {
// default. simple intersection only
if (_writeA == false && _writeB == false && _writeOverlap == false) {
_bedA->reportBedRangeNewLine(a,s,e);
}
// -wa -wbwrite the original A and B
else if (_writeA == true && _writeB == true) {
_bedA->reportBedTab(a);
_bedB->reportBedNewLine(b);
}
// -wa write just the original A
else if (_writeA == true) {
_bedA->reportBedNewLine(a);
}
// -wb write the intersected portion of A and the original B
else if (_writeB == true) {
_bedA->reportBedRangeTab(a,s,e);
_bedB->reportBedNewLine(b);
}
// -wo write the original A and B plus the no. of overlapping bases.
else if (_writeOverlap == true) {
_bedA->reportBedTab(a);
_bedB->reportBedTab(b);
printf("%d\n", overlapBases);
}
}
void BedIntersect::ReportOverlapSummary(const BED &a, const int &numOverlapsFound) {
// -u just report the fact that there was >= 1 overlaps
if (_anyHit && (numOverlapsFound >= 1)) {
_bedA->reportBedNewLine(a);
}
// -c report the total number of features overlapped in B
else if (_writeCount) {
_bedA->reportBedTab(a);
printf("%d\n", numOverlapsFound);
}
// -v report iff there were no overlaps
else if (_noHit && (numOverlapsFound == 0)) {
_bedA->reportBedNewLine(a);
}
// -wao the user wants to force the reporting of 0 overlap
else if (_writeAllOverlap && (numOverlapsFound == 0)) {
_bedA->reportBedTab(a);
_bedB->reportNullBedTab();
printf("0\n");
}
}
void BedIntersect::Scan(BED *x, vector<BED *> *windowX, BedLineStatus xStatus,
const BED &y, vector<BED *> *windowY, BedLineStatus yStatus) {
if (xStatus != BED_VALID) {
return;
}
std::vector<BED *>::iterator wYIter = windowY->begin();
while (wYIter != windowY->end()) {
if (leftOf(*(*wYIter), *x) == true) {
(*wYIter)->finished = true;
wYIter = windowY->erase(wYIter); // erase auto-increments to the next position
}
else if (overlaps((*wYIter)->start, (*wYIter)->end, x->start, x->end) > 0) {
if (_lastPick == 0) {
AddHits(x, *(*wYIter));
}
else {
AddHits(*wYIter, *x);
}
++wYIter; // force incrementing
}
}
if (leftOf(*x,y) == false)
windowX->push_back(x);
else {
x->finished = true;
}
// dump the buffered results (if any)
FlushOutputBuffer();
}
void BedIntersect::AddHits(BED *x, const BED &y) {
if (x->overlaps.empty() == true)
_outputBuffer.push(x);
x->overlaps.push_back(y);
}
void BedIntersect::FlushOutputBuffer(bool final) {
while (_outputBuffer.empty() == false)
{
if (final == false && _outputBuffer.front()->finished == false)
break;
processHits(*_outputBuffer.front(), _outputBuffer.front()->overlaps);
// remove the finished BED entry from the heap
delete _outputBuffer.front();
_outputBuffer.pop();
}
}
vector<BED*>* BedIntersect::GetWindow(const string &chrom, bool isA) {
// iterator to test if a window for a given chrom exists.
map<string, vector<BED*> >::iterator it;
// grab the current window for A or B, depending on
// the request. if a window hasn't yet been created
// for the requested chrom, create one.
if (isA) {
it = _windowA.find(chrom);
if (it != _windowA.end()) {
return & _windowA[chrom];
}
else {
_windowA.insert(pair<string, vector<BED *> >(chrom, vector<BED *>()));
return & _windowA[chrom];
}
}
else {
it = _windowB.find(chrom);
if (it != _windowB.end()) {
return & _windowB[chrom];
}
else {
_windowB.insert(pair<string, vector<BED *> >(chrom, vector<BED *>()));
return & _windowB[chrom];
}
}
}
void BedIntersect::ChromSwitch(const string &chrom) {
vector<BED*>::iterator windowAIter = _windowA[chrom].begin();
vector<BED*>::iterator windowAEnd = _windowA[chrom].end();
for (; windowAIter != windowAEnd; ++windowAIter)
(*windowAIter)->finished = true;
vector<BED*>::iterator windowBIter = _windowB[chrom].begin();
vector<BED*>::iterator windowBEnd = _windowB[chrom].end();
for (; windowBIter != windowBEnd; ++windowBIter)
(*windowBIter)->finished = true;
FlushOutputBuffer();
}
void BedIntersect::IntersectBed() {
int aLineNum = 0;
int bLineNum = 0;
// current feature from each file
BED *a, *b, *prevA, *prevB;
// status of the current lines
BedLineStatus aStatus, bStatus;
// open the files; get the first line from each
_bedA->Open();
_bedB->Open();
prevA = NULL;
prevB = NULL;
a = new BED();
b = new BED();
aStatus = _bedA->GetNextBed(*a, aLineNum);
bStatus = _bedB->GetNextBed(*b, bLineNum);
cout << a->chrom << " " << a->start << " " << a->chrom << " " << b->start << endl;
while (aStatus != BED_INVALID || bStatus != BED_INVALID) {
if ((a->start <= b->start) && (a->chrom == b->chrom)) {
prevA = a;
_lastPick = 0;
Scan(a, GetWindow(a->chrom, true), aStatus,
*b, GetWindow(a->chrom, false), bStatus);
a = new BED();
aStatus = _bedA->GetNextBed(*a, aLineNum);
}
else if ((a->start > b->start) && (a->chrom == b->chrom)) {
prevB = b;
_lastPick = 1;
Scan(b, GetWindow(b->chrom, false), bStatus,
*a, GetWindow(b->chrom, true), aStatus);
b = new BED();
bStatus = _bedB->GetNextBed(*b, bLineNum);
}
else if (a->chrom != b->chrom) {
// A was most recently read
if (_lastPick == 0) {
prevB = b;
while (b->chrom == prevA->chrom){
_windowB[prevA->chrom].push_back(b);
b = new BED();
bStatus = _bedB->GetNextBed(*b, bLineNum);
}
Scan(prevA, GetWindow(prevA->chrom, true), aStatus,
*prevB, GetWindow(prevA->chrom, false), bStatus);
}
// B was most recently read
else {
prevA = a;
while (a->chrom == prevB->chrom) {
_windowA[prevB->chrom].push_back(a);
a = new BED();
aStatus = _bedA->GetNextBed(*a, aLineNum);
}
Scan(prevB, GetWindow(prevB->chrom, false), bStatus,
*prevA, GetWindow(prevB->chrom, true), aStatus);
}
FlushOutputBuffer(true);
}
if (prevA!=NULL&&prevB!=NULL)
//cout << prevA->chrom << " " << a->chrom << " " << a->start << " "
// << prevB->chrom << " " << b->chrom << " " << b->start << "\n";
if (aStatus == BED_INVALID) a->start = INT_MAX;
if (bStatus == BED_INVALID) b->start = INT_MAX;
}
// clear out the final bit of staged output
FlushOutputBuffer(true);
// close the files
_bedA->Close();
_bedB->Close();
}
/*****************************************************************************
intersectBed.h
(c) 2009 - Aaron Quinlan
Hall Laboratory
Department of Biochemistry and Molecular Genetics
University of Virginia
aaronquinlan@gmail.com
Licenced under the GNU General Public License 2.0 license.
******************************************************************************/
#ifndef INTERSECTBED_H
#define INTERSECTBED_H
#include "bedFile.h"
// #include "BamReader.h"
// #include "BamWriter.h"
// #include "BamAncillary.h"
// #include "BamAux.h"
// using namespace BamTools;
#include <vector>
#include <queue>
#include <iostream>
#include <fstream>
#include <stdlib.h>
using namespace std;
class BedIntersect {
public:
// constructor
BedIntersect(string bedAFile, string bedBFile, bool anyHit,
bool writeA, bool writeB, bool writeOverlap, bool writeAllOverlap,
float overlapFraction, bool noHit, bool writeCount, bool forceStrand,
bool reciprocal, bool obeySplits, bool bamInput, bool bamOutput);
// destructor
~BedIntersect(void);
private:
//------------------------------------------------
// private attributes
//------------------------------------------------
string _bedAFile;
string _bedBFile;
bool _writeA; // should the original A feature be reported?
bool _writeB; // should the original B feature be reported?
bool _writeOverlap;
bool _writeAllOverlap;
bool _forceStrand;
bool _reciprocal;
float _overlapFraction;
bool _anyHit;
bool _noHit;
bool _writeCount; // do we want a count of the number of overlaps in B?
bool _obeySplits;
bool _bamInput;
bool _bamOutput;
bool _printable;
queue<BED*> _outputBuffer;
bool _lastPick;
map<string, vector<BED*> > _windowA;
map<string, vector<BED*> > _windowB;
// instance of a bed file class.
BedFile *_bedA, *_bedB;
//------------------------------------------------
// private methods
//------------------------------------------------
void IntersectBed(istream &bedInput);
void Scan(BED *x, vector<BED *> *windowX, BedLineStatus xStatus,
const BED &y, vector<BED *> *windowY, BedLineStatus yStatus);
void AddHits(BED *x, const BED &y);
void FlushOutputBuffer(bool final = false);
vector<BED*>* GetWindow(const string &chrom, bool isA);
void ChromSwitch(const string &chrom);
void IntersectBed();
void IntersectBam(string bamFile);
bool processHits(BED &a, vector<BED> &hits);
bool FindOverlaps(const BED &a, vector<BED> &hits);
bool FindOneOrMoreOverlap(const BED &a);
void ReportOverlapDetail(const int &overlapBases, const BED &a, const BED &b,
const CHRPOS &s, const CHRPOS &e);
void ReportOverlapSummary(const BED &a, const int &numOverlapsFound);
void ReportHits(set<BED> &A, set<BED> &B);
};
#endif /* INTERSECTBED_H */
/*****************************************************************************
intersectMain.cpp
(c) 2009 - Aaron Quinlan
Hall Laboratory
Department of Biochemistry and Molecular Genetics
University of Virginia
aaronquinlan@gmail.com
Licenced under the GNU General Public License 2.0 license.
******************************************************************************/
#include "fjoin.h"
#include "version.h"
using namespace std;
// define our program name
#define PROGRAM_NAME "fjoin"
// define our parameter checking macro
#define PARAMETER_CHECK(param, paramLen, actualLen) (strncmp(argv[i], param, min(actualLen, paramLen))== 0) && (actualLen == paramLen)
// function declarations
void ShowHelp(void);
int main(int argc, char* argv[]) {
// our configuration variables
bool showHelp = false;
// input files
string bedAFile;
string bedBFile;
// input arguments
float overlapFraction = 1E-9;
bool haveBedA = false;
bool haveBedB = false;
bool noHit = false;
bool anyHit = false;
bool writeA = false;
bool writeB = false;
bool writeCount = false;
bool writeOverlap = false;
bool writeAllOverlap = false;
bool haveFraction = false;
bool reciprocalFraction = false;
bool forceStrand = false;
bool obeySplits = false;
bool inputIsBam = false;
bool outputIsBam = true;
// check to see if we should print out some help
if(argc <= 1) showHelp = true;
for(int i = 1; i < argc; i++) {
int parameterLength = (int)strlen(argv[i]);
if((PARAMETER_CHECK("-h", 2, parameterLength)) ||
(PARAMETER_CHECK("--help", 5, parameterLength))) {
showHelp = true;
}
}
if(showHelp) ShowHelp();
// do some parsing (all of these parameters require 2 strings)
for(int i = 1; i < argc; i++) {
int parameterLength = (int)strlen(argv[i]);
if(PARAMETER_CHECK("-a", 2, parameterLength)) {
if ((i+1) < argc) {
haveBedA = true;
outputIsBam = false;
bedAFile = argv[i + 1];
i++;
}
}
else if(PARAMETER_CHECK("-abam", 5, parameterLength)) {
if ((i+1) < argc) {
haveBedA = true;
inputIsBam = true;
bedAFile = argv[i + 1];
i++;
}
}
else if(PARAMETER_CHECK("-b", 2, parameterLength)) {
if ((i+1) < argc) {
haveBedB = true;
bedBFile = argv[i + 1];
i++;
}
}
else if(PARAMETER_CHECK("-bed", 4, parameterLength)) {
outputIsBam = false;
}
else if(PARAMETER_CHECK("-u", 2, parameterLength)) {
anyHit = true;
}
else if(PARAMETER_CHECK("-f", 2, parameterLength)) {
if ((i+1) < argc) {
haveFraction = true;
overlapFraction = atof(argv[i + 1]);
i++;
}
}
else if(PARAMETER_CHECK("-wa", 3, parameterLength)) {
writeA = true;
}
else if(PARAMETER_CHECK("-wb", 3, parameterLength)) {
writeB = true;
}
else if(PARAMETER_CHECK("-wo", 3, parameterLength)) {
writeOverlap = true;
}
else if(PARAMETER_CHECK("-wao", 4, parameterLength)) {
writeAllOverlap = true;
writeOverlap = true;
}
else if(PARAMETER_CHECK("-c", 2, parameterLength)) {
writeCount = true;
}
else if(PARAMETER_CHECK("-r", 2, parameterLength)) {
reciprocalFraction = true;
}
else if (PARAMETER_CHECK("-v", 2, parameterLength)) {
noHit = true;
}
else if (PARAMETER_CHECK("-s", 2, parameterLength)) {
forceStrand = true;
}
else if (PARAMETER_CHECK("-split", 6, parameterLength)) {
obeySplits = true;
}
else {
cerr << endl << "*****ERROR: Unrecognized parameter: " << argv[i] << " *****" << endl << endl;
showHelp = true;
}
}
// make sure we have both input files
if (!haveBedA || !haveBedB) {
cerr << endl << "*****" << endl << "*****ERROR: Need -a and -b files. " << endl << "*****" << endl;
showHelp = true;
}
if (anyHit && noHit) {
cerr << endl << "*****" << endl << "*****ERROR: Request either -u OR -v, not both." << endl << "*****" << endl;
showHelp = true;
}
if (writeB && writeCount) {
cerr << endl << "*****" << endl << "*****ERROR: Request either -wb OR -c, not both." << endl << "*****" << endl;
showHelp = true;
}
if (writeCount && writeOverlap) {
cerr << endl << "*****" << endl << "*****ERROR: Request either -wb OR -wo, not both." << endl << "*****" << endl;
showHelp = true;
}
if (writeA && writeOverlap) {
cerr << endl << "*****" << endl << "*****ERROR: Request either -wa OR -wo, not both." << endl << "*****" << endl;
showHelp = true;
}
if (writeB && writeOverlap) {
cerr << endl << "*****" << endl << "*****ERROR: Request either -wb OR -wo, not both." << endl << "*****" << endl;
showHelp = true;
}
if (reciprocalFraction && !haveFraction) {
cerr << endl << "*****" << endl << "*****ERROR: If using -r, you need to define -f." << endl << "*****" << endl;
showHelp = true;
}
if (anyHit && writeCount) {
cerr << endl << "*****" << endl << "*****ERROR: Request either -u OR -c, not both." << endl << "*****" << endl;
showHelp = true;
}
if (anyHit && writeB) {
cerr << endl << "*****" << endl << "*****ERROR: Request either -u OR -wb, not both." << endl << "*****" << endl;
showHelp = true;
}
if (anyHit && writeOverlap) {
cerr << endl << "*****" << endl << "*****ERROR: Request either -u OR -wo, not both." << endl << "*****" << endl;
showHelp = true;
}
if (!showHelp) {
BedIntersect *bi = new BedIntersect(bedAFile, bedBFile, anyHit, writeA, writeB, writeOverlap,
writeAllOverlap, overlapFraction, noHit, writeCount, forceStrand,
reciprocalFraction, obeySplits, inputIsBam, outputIsBam);
delete bi;
return 0;
}
else {
ShowHelp();
}
}
void ShowHelp(void) {
cerr << endl << "Program: " << PROGRAM_NAME << " (v" << VERSION << ")" << endl;
cerr << "Author: Aaron Quinlan (aaronquinlan@gmail.com)" << endl;
cerr << "Summary: Report overlaps between two feature files." << endl << endl;
cerr << "Usage: " << PROGRAM_NAME << " [OPTIONS] -a <bed/gff/vcf> -b <bed/gff/vcf>" << endl << endl;
cerr << "Options: " << endl;
cerr << "\t-abam\t" << "The A input file is in BAM format. Output will be BAM as well." << endl << endl;
cerr << "\t-bed\t" << "When using BAM input (-abam), write output as BED. The default" << endl;
cerr << "\t\tis to write output in BAM when using -abam." << endl << endl;
cerr << "\t-wa\t" << "Write the original entry in A for each overlap." << endl << endl;
cerr << "\t-wb\t" << "Write the original entry in B for each overlap." << endl;
cerr << "\t\t- Useful for knowing _what_ A overlaps. Restricted by -f and -r." << endl << endl;
cerr << "\t-wo\t" << "Write the original A and B entries plus the number of base" << endl;
cerr << "\t\tpairs of overlap between the two features." << endl;
cerr << "\t\t- Overlaps restricted by -f and -r." << endl;
cerr << "\t\t Only A features with overlap are reported." << endl << endl;
cerr << "\t-wao\t" << "Write the original A and B entries plus the number of base" << endl;
cerr << "\t\tpairs of overlap between the two features." << endl;
cerr << "\t\t- Overlapping features restricted by -f and -r." << endl;
cerr << "\t\t However, A features w/o overlap are also reported" << endl;
cerr << "\t\t with a NULL B feature and overlap = 0." << endl << endl;
cerr << "\t-u\t" << "Write the original A entry _once_ if _any_ overlaps found in B." << endl;
cerr << "\t\t- In other words, just report the fact >=1 hit was found." << endl;
cerr << "\t\t- Overlaps restricted by -f and -r." << endl << endl;
cerr << "\t-c\t" << "For each entry in A, report the number of overlaps with B." << endl;
cerr << "\t\t- Reports 0 for A entries that have no overlap with B." << endl;
cerr << "\t\t- Overlaps restricted by -f and -r." << endl << endl;
cerr << "\t-v\t" << "Only report those entries in A that have _no overlaps_ with B." << endl;
cerr << "\t\t- Similar to \"grep -v\" (an homage)." << endl << endl;
cerr << "\t-f\t" << "Minimum overlap required as a fraction of A." << endl;
cerr << "\t\t- Default is 1E-9 (i.e., 1bp)." << endl;
cerr << "\t\t- FLOAT (e.g. 0.50)" << endl << endl;
cerr << "\t-r\t" << "Require that the fraction overlap be reciprocal for A and B." << endl;
cerr << "\t\t- In other words, if -f is 0.90 and -r is used, this requires" << endl;
cerr << "\t\t that B overlap 90% of A and A _also_ overlaps 90% of B." << endl << endl;
cerr << "\t-s\t" << "Force strandedness. That is, only report hits in B that" << endl;
cerr << "\t\toverlap A on the same strand." << endl;
cerr << "\t\t- By default, overlaps are reported without respect to strand." << endl << endl;
cerr << "\t-split\t" << "Treat \"split\" BAM or BED12 entries as distinct BED intervals." << endl << endl;
// end the program here
exit(1);
}
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