From 520ccaaef27d06bd6d2a5575592687c10e57f27e Mon Sep 17 00:00:00 2001
From: Aaron <aaronquinlan@gmail.com>
Date: Wed, 7 Dec 2011 21:31:00 -0500
Subject: [PATCH] Use new GetNextBed() interface.
---
src/cuffToTrans/Makefile | 44 -----
src/fjoin/Makefile | 42 -----
src/fjoin/fjoin.cpp | 350 ---------------------------------------
src/fjoin/fjoin.h | 114 -------------
src/fjoin/fjoinMain.cpp | 271 ------------------------------
5 files changed, 821 deletions(-)
delete mode 100644 src/cuffToTrans/Makefile
delete mode 100644 src/fjoin/Makefile
delete mode 100644 src/fjoin/fjoin.cpp
delete mode 100644 src/fjoin/fjoin.h
delete mode 100644 src/fjoin/fjoinMain.cpp
diff --git a/src/cuffToTrans/Makefile b/src/cuffToTrans/Makefile
deleted file mode 100644
index 6e417369..00000000
--- a/src/cuffToTrans/Makefile
+++ /dev/null
@@ -1,44 +0,0 @@
-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
diff --git a/src/fjoin/Makefile b/src/fjoin/Makefile
deleted file mode 100644
index c98c6696..00000000
--- a/src/fjoin/Makefile
+++ /dev/null
@@ -1,42 +0,0 @@
-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
diff --git a/src/fjoin/fjoin.cpp b/src/fjoin/fjoin.cpp
deleted file mode 100644
index b1de416a..00000000
--- a/src/fjoin/fjoin.cpp
+++ /dev/null
@@ -1,350 +0,0 @@
-/*****************************************************************************
- 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();
-}
-
-
diff --git a/src/fjoin/fjoin.h b/src/fjoin/fjoin.h
deleted file mode 100644
index c7aabd46..00000000
--- a/src/fjoin/fjoin.h
+++ /dev/null
@@ -1,114 +0,0 @@
-/*****************************************************************************
- 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 */
diff --git a/src/fjoin/fjoinMain.cpp b/src/fjoin/fjoinMain.cpp
deleted file mode 100644
index 48a2ad98..00000000
--- a/src/fjoin/fjoinMain.cpp
+++ /dev/null
@@ -1,271 +0,0 @@
-/*****************************************************************************
- 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);
-
-}
--
GitLab