diff --git a/Makefile b/Makefile
index d61800288ed348916f91fa143c665f5303bd7b06..a162fbd61378583c83dfd437dae47ce5b7a7e87e 100755
--- a/Makefile
+++ b/Makefile
@@ -19,7 +19,7 @@ export CXX ?= g++
#include includes/$(BLD_PLATFORM).inc
# define our source subdirectories
-SUBDIRS = $(SRC_DIR)/closestBed $(SRC_DIR)/complementBed $(SRC_DIR)/coverageBed $(SRC_DIR)/intersectBed $(SRC_DIR)/mergeBed $(SRC_DIR)/genomeCoverageBed $(SRC_DIR)/fastaFromBed $(SRC_DIR)/sortBed $(SRC_DIR)/windowBed $(SRC_DIR)/subtractBed $(SRC_DIR)/linksBed $(SRC_DIR)/peIntersectBed
+SUBDIRS = $(SRC_DIR)/closestBed $(SRC_DIR)/complementBed $(SRC_DIR)/coverageBed $(SRC_DIR)/intersectBed $(SRC_DIR)/mergeBed $(SRC_DIR)/genomeCoverageBed $(SRC_DIR)/fastaFromBed $(SRC_DIR)/shuffleBed $(SRC_DIR)/slopBed $(SRC_DIR)/sortBed $(SRC_DIR)/windowBed $(SRC_DIR)/subtractBed $(SRC_DIR)/linksBed $(SRC_DIR)/peIntersectBed $(SRC_DIR)/pairToPair $(SRC_DIR)/maskFastaFromBed
UTIL_SUBDIRS = $(SRC_DIR)/utils/lineFileUtilities $(SRC_DIR)/utils/bedFile $(SRC_DIR)/utils/bedFilePE $(SRC_DIR)/utils/sequenceUtilities
diff --git a/bin/closestBed b/bin/closestBed
index 1ca82d27dec77e6907bca3e05c76489636e22e00..07be060e28c2894d27dc19f2aadca549c328af70 100755
Binary files a/bin/closestBed and b/bin/closestBed differ
diff --git a/bin/complementBed b/bin/complementBed
index 04300802881a7b9c016e2393a982a287316fc053..9a296e6ce7529526d89ad1de7c74a280588438af 100755
Binary files a/bin/complementBed and b/bin/complementBed differ
diff --git a/bin/coverageBed b/bin/coverageBed
index edf902cd4d4287284de7afe7f84cd50b7aead0d2..2f39bbb0b1a537d91324ee34c76cd47f40d09f16 100755
Binary files a/bin/coverageBed and b/bin/coverageBed differ
diff --git a/bin/fastaFromBed b/bin/fastaFromBed
index 1e81d12558f1a586cad7afa91fbe50df25f352bd..18ffbe080cd8011745e9037bc8008c1711830b9c 100755
Binary files a/bin/fastaFromBed and b/bin/fastaFromBed differ
diff --git a/bin/genomeCoverageBed b/bin/genomeCoverageBed
index 5c3db43246254d0e507f5a2c13d54ec09a78e416..9a43ab83188f25aede5d981ed70e7e6ff55f7355 100755
Binary files a/bin/genomeCoverageBed and b/bin/genomeCoverageBed differ
diff --git a/bin/linksBed b/bin/linksBed
index 26d33cad1dbda6d31379ec32536903742e5215c6..b501ef8bd490d1c3a8d3b27be2ffd953ed238f8b 100755
Binary files a/bin/linksBed and b/bin/linksBed differ
diff --git a/bin/maskFastaFromBed b/bin/maskFastaFromBed
new file mode 100755
index 0000000000000000000000000000000000000000..f1c9b08b09b61b5cc31a7a0fd22049ae0ab12894
Binary files /dev/null and b/bin/maskFastaFromBed differ
diff --git a/bin/mergeBed b/bin/mergeBed
index 12d3396e0115962768cfeda6f7e1c9c90b6dbbeb..c6f54ae46e8340ac66a8bbacfede5c18171d6e2e 100755
Binary files a/bin/mergeBed and b/bin/mergeBed differ
diff --git a/bin/pairToPair b/bin/pairToPair
new file mode 100755
index 0000000000000000000000000000000000000000..1f31950a3f6e6b007e56de95fc3476adf6fce092
Binary files /dev/null and b/bin/pairToPair differ
diff --git a/bin/peIntersectBed b/bin/peIntersectBed
index c006038cfd166d23b8a63335dbde94ccd38916c0..4c39e977a93b05ffe7d369d924db6e26bc93159e 100755
Binary files a/bin/peIntersectBed and b/bin/peIntersectBed differ
diff --git a/bin/shuffleBed b/bin/shuffleBed
new file mode 100755
index 0000000000000000000000000000000000000000..68d55d71049856aa861c4e49151d742418cbc3e3
Binary files /dev/null and b/bin/shuffleBed differ
diff --git a/bin/slopBed b/bin/slopBed
new file mode 100755
index 0000000000000000000000000000000000000000..c2c7c3bee4048ca5e045d5400819817752adc7b2
Binary files /dev/null and b/bin/slopBed differ
diff --git a/bin/sortBed b/bin/sortBed
index 9b6ee50d4d040d05d132c4b36aa65ce1c86211fd..71379afb5216dcbe595f02b70aa190260ea5919c 100755
Binary files a/bin/sortBed and b/bin/sortBed differ
diff --git a/bin/subtractBed b/bin/subtractBed
index 7014f75affd723cadaf17ee4434792a2f5daac0d..9ad38b265e5f8d6004f9c1e6ca93936b843f8df6 100755
Binary files a/bin/subtractBed and b/bin/subtractBed differ
diff --git a/bin/windowBed b/bin/windowBed
index fa1de756abb575301308cbd0ac8ab59fde408a77..075e1b68c2cb95e81daef6c66b38b20509daf5d0 100755
Binary files a/bin/windowBed and b/bin/windowBed differ
diff --git a/obj/bedFilePE.o b/obj/bedFilePE.o
index 6b65a11ad79eb2758ba64887cd0021fdd8519549..66c0da3548b8a6b566d91e135ecbd639b57342a5 100644
Binary files a/obj/bedFilePE.o and b/obj/bedFilePE.o differ
diff --git a/obj/closestBed.o b/obj/closestBed.o
index dfbf99a33f0211218030550cfa6af7b3e648767b..0af6ce2e0d70fe5e7993a59fccb37107ff4e9ad9 100644
Binary files a/obj/closestBed.o and b/obj/closestBed.o differ
diff --git a/obj/closestMain.o b/obj/closestMain.o
index d912aced832eda6140e8fd221f25970710ccb57e..ce51ad95642bc9d350677e186d23401e3c5a9827 100644
Binary files a/obj/closestMain.o and b/obj/closestMain.o differ
diff --git a/obj/complementMain.o b/obj/complementMain.o
index 1c8ce5c50bc8f83e04b891b6d512949f7eaa307f..2c6dfe9ca3dca66253e68b5aa9a29a3929d194fc 100644
Binary files a/obj/complementMain.o and b/obj/complementMain.o differ
diff --git a/obj/coverageBed.o b/obj/coverageBed.o
index 5f5d66a692e2040c0af50d3ffc7e05dfdab7ca4c..f120e6f002aba7772f51358860985f13ab61f2a2 100644
Binary files a/obj/coverageBed.o and b/obj/coverageBed.o differ
diff --git a/obj/coverageMain.o b/obj/coverageMain.o
index 607976d8eaba31e7c3afd0047a82d9bc415f2bbc..5f4963c94c38734237b4c0fa6091a7d773a323d6 100644
Binary files a/obj/coverageMain.o and b/obj/coverageMain.o differ
diff --git a/obj/fastaFromBed.o b/obj/fastaFromBed.o
index be08d863a8d91062baf31455122dd764815ccf71..260de5fa368aef290a5b49d174fd68c7bf48143f 100644
Binary files a/obj/fastaFromBed.o and b/obj/fastaFromBed.o differ
diff --git a/obj/fastaFromBedMain.o b/obj/fastaFromBedMain.o
index d5d645ff1fa614ad434b8f1d42fad852e110a396..fd464f4e9107131cd18a0e13cdefb9b973e77049 100644
Binary files a/obj/fastaFromBedMain.o and b/obj/fastaFromBedMain.o differ
diff --git a/obj/genomeCoverageBed.o b/obj/genomeCoverageBed.o
index 9b27e8d025835d6eceab0f9709cb5994690bb6c1..796ab559e8cb21fb596350f2e7420b306fb75bd2 100644
Binary files a/obj/genomeCoverageBed.o and b/obj/genomeCoverageBed.o differ
diff --git a/obj/genomeCoverageMain.o b/obj/genomeCoverageMain.o
index 2ab301eb359a396fa6f7915358c9f26c53ec4128..901de722a41e4fa24eb2429bf80c205ba93d52f5 100644
Binary files a/obj/genomeCoverageMain.o and b/obj/genomeCoverageMain.o differ
diff --git a/obj/linksBed.o b/obj/linksBed.o
index ec65f8caab724e44b2a4e20676bc43ef3328b6b7..a8a7f07e35b91a05d9070317a7b8be242fd1bb69 100644
Binary files a/obj/linksBed.o and b/obj/linksBed.o differ
diff --git a/obj/linksMain.o b/obj/linksMain.o
index e9f1d97599557b7f6fb6ec06f0354a34bfe0f943..4c5f8c538fc10308f909b3d10b3da5a46bb36e7e 100644
Binary files a/obj/linksMain.o and b/obj/linksMain.o differ
diff --git a/obj/maskFastaFromBed.o b/obj/maskFastaFromBed.o
new file mode 100644
index 0000000000000000000000000000000000000000..adb815cf587e449cb882d1440c5d5bb17dc0c63b
Binary files /dev/null and b/obj/maskFastaFromBed.o differ
diff --git a/obj/maskFastaFromBedMain.o b/obj/maskFastaFromBedMain.o
new file mode 100644
index 0000000000000000000000000000000000000000..8b50539aa1f7b67f4aa7104d6de3b10abc7810b0
Binary files /dev/null and b/obj/maskFastaFromBedMain.o differ
diff --git a/obj/pairToPair.o b/obj/pairToPair.o
new file mode 100644
index 0000000000000000000000000000000000000000..2cfb188d7c9590758ba956b3a88edce9e9fa4c26
Binary files /dev/null and b/obj/pairToPair.o differ
diff --git a/obj/pairToPairMain.o b/obj/pairToPairMain.o
new file mode 100644
index 0000000000000000000000000000000000000000..ae50cccada0338825712c4b9ffbeaa839f329336
Binary files /dev/null and b/obj/pairToPairMain.o differ
diff --git a/obj/peIntersectBed.o b/obj/peIntersectBed.o
index fe06952a59bfb48e5fdb32b75f10cffb47ecb810..c95bda37a96b505ade3a83cfe3f49bec12fd0425 100644
Binary files a/obj/peIntersectBed.o and b/obj/peIntersectBed.o differ
diff --git a/obj/peIntersectMain.o b/obj/peIntersectMain.o
index 27fa3b4b6b4db6bae007132b87141050393cc89e..92de575c95d843864c13e75069a35f4101f4f35d 100644
Binary files a/obj/peIntersectMain.o and b/obj/peIntersectMain.o differ
diff --git a/obj/sequenceUtils.o b/obj/sequenceUtils.o
index ba80be57f249033b0641e57abd309e12bf316c29..38fe6e47944f70a7c672162231bb8c0b54814910 100644
Binary files a/obj/sequenceUtils.o and b/obj/sequenceUtils.o differ
diff --git a/obj/shuffleBed.o b/obj/shuffleBed.o
new file mode 100644
index 0000000000000000000000000000000000000000..e1852b36a38d1e4fda2701efad8e232f15e98193
Binary files /dev/null and b/obj/shuffleBed.o differ
diff --git a/obj/shuffleBedMain.o b/obj/shuffleBedMain.o
new file mode 100644
index 0000000000000000000000000000000000000000..0c327e88b80654efab061ccc7cca10cba0b897fb
Binary files /dev/null and b/obj/shuffleBedMain.o differ
diff --git a/obj/slopBed.o b/obj/slopBed.o
new file mode 100644
index 0000000000000000000000000000000000000000..e04730d841d9213809b81d308bd95828bb0aa391
Binary files /dev/null and b/obj/slopBed.o differ
diff --git a/obj/slopBedMain.o b/obj/slopBedMain.o
new file mode 100644
index 0000000000000000000000000000000000000000..ebf824fa030386ee1c2e2bc6164900363b67d835
Binary files /dev/null and b/obj/slopBedMain.o differ
diff --git a/obj/sortBed.o b/obj/sortBed.o
index 5686d0e4964d5f7d039ee21971c431bfe4e8869a..22b7626f547b3e6f79e656829d1d7d566a720a46 100644
Binary files a/obj/sortBed.o and b/obj/sortBed.o differ
diff --git a/obj/sortMain.o b/obj/sortMain.o
index 3c8d8f56d10bf30caace87f203a30481e72922f2..ec389c8d9d4b610bfd78a4a2eaecbdac1199f285 100644
Binary files a/obj/sortMain.o and b/obj/sortMain.o differ
diff --git a/obj/subtractBed.o b/obj/subtractBed.o
index bf06bf9a21591b9514b0219ffa57445c66c532b9..bc8768e677bfe8e74a2c36a77441c9383597bd1d 100644
Binary files a/obj/subtractBed.o and b/obj/subtractBed.o differ
diff --git a/obj/subtractMain.o b/obj/subtractMain.o
index 577a06362ec7db267a43f236f02e2be4ec069611..e2be0d73e8851e959fec5a5787e03ba78fe0e33f 100644
Binary files a/obj/subtractMain.o and b/obj/subtractMain.o differ
diff --git a/obj/windowBed.o b/obj/windowBed.o
index ba94fd6e689bf811c2ed366c252a76735188f38f..5c043aadd66b3116b996d6b931d282cbe1d956de 100644
Binary files a/obj/windowBed.o and b/obj/windowBed.o differ
diff --git a/obj/windowMain.o b/obj/windowMain.o
index 7e513d449ea12ea5a6bcebcd053f87af9a0bc65f..0ae2e3b64278c545802121d3402f411f537dd9e2 100644
Binary files a/obj/windowMain.o and b/obj/windowMain.o differ
diff --git a/src/pairToPair/Makefile b/src/pairToPair/Makefile
new file mode 100755
index 0000000000000000000000000000000000000000..248c635c88aca5a03adb3a0db999a78df3e13ba5
--- /dev/null
+++ b/src/pairToPair/Makefile
@@ -0,0 +1,46 @@
+CXX = g++
+CXXFLAGS = -O3 -Wall
+LDFLAGS =
+
+UTILITIES_DIR = ../utils/
+OBJ_DIR = ../../obj/
+BIN_DIR = ../../bin/
+
+# -------------------
+# define our includes
+# -------------------
+INCLUDES = -I$(UTILITIES_DIR)/bedFilePE/ -I$(UTILITIES_DIR)/bedFile/ -I$(UTILITIES_DIR)/lineFileUtilities/ -I$(UTILITIES_DIR)/version/
+
+# ----------------------------------
+# define our source and object files
+# ----------------------------------
+SOURCES= pairToPairMain.cpp pairToPair.cpp
+OBJECTS= $(SOURCES:.cpp=.o)
+_EXT_OBJECTS=bedFilePE.o bedFile.o lineFileUtilities.o
+EXT_OBJECTS=$(patsubst %,$(OBJ_DIR)/%,$(_EXT_OBJECTS))
+BUILT_OBJECTS= $(patsubst %,$(OBJ_DIR)/%,$(OBJECTS))
+PROGRAM= pairToPair
+LIBS=
+
+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)/bedFilePE/
+ @$(MAKE) --no-print-directory -C $(UTILITIES_DIR)/lineFileUtilities/
+
+clean:
+ @echo "Cleaning up."
+ @rm -f $(OBJ_DIR)/* $(BIN_DIR)/*
+
+.PHONY: clean
diff --git a/src/pairToPair/pairToPair.cpp b/src/pairToPair/pairToPair.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..118fd3fdcb3bef0ef7e67f8ae88f4cef93780bcd
--- /dev/null
+++ b/src/pairToPair/pairToPair.cpp
@@ -0,0 +1,226 @@
+/*
+ BEDTools: pairToPair.cpp
+
+ Created by Aaron Quinlan Spring 2009.
+ Copyright 2009 Aaron Quinlan. All rights reserved.
+
+ Summary: Looks for overlaps between paired-end reads / SVs (BEDPE format) and a BED file.
+*/
+#include "lineFileUtilities.h"
+#include "pairToPair.h"
+
+
+/*
+ Constructor
+*/
+PairToPair::PairToPair(string &bedAFilePE, string &bedBFilePE, float &overlapFraction,
+ string &searchType) {
+
+ this->bedAFilePE = bedAFilePE;
+ this->bedBFilePE = bedBFilePE;
+ this->overlapFraction = overlapFraction;
+ this->searchType = searchType;
+
+ this->bedA = new BedFilePE(bedAFilePE);
+ this->bedB = new BedFilePE(bedBFilePE);
+}
+
+
+/*
+ Destructor
+*/
+PairToPair::~PairToPair(void) {
+}
+
+
+
+void PairToPair::IntersectPairs() {
+
+ // load the "B" bed file into a map so
+ // that we can easily compare "A" to it for overlaps
+ bedB->loadBedPEFileIntoMap();
+
+ string bedLine;
+ int lineNum = 0;
+
+ // are we dealing with a file?
+ if (bedA->bedFile != "stdin") {
+
+ // open the BEDPE file for reading
+ ifstream bed(bedA->bedFile.c_str(), ios::in);
+ if ( !bed ) {
+ cerr << "Error: The requested bedpe file (" <<bedA->bedFile << ") could not be opened. Exiting!" << endl;
+ exit (1);
+ }
+
+ BEDPE a;
+ // process each entry in A
+ while (getline(bed, bedLine)) {
+
+ // split the current line into ditinct fields
+ vector<string> bedFields;
+ Tokenize(bedLine,bedFields);
+
+ lineNum++;
+
+ // find the overlaps with B if it's a valid BED entry.
+ if (bedA->parseBedPELine(a, bedFields, lineNum)) {
+ vector<BED> hitsA1B1, hitsA1B2, hitsA2B1, hitsA2B2;
+
+ FindOverlaps(a, hitsA1B1, hitsA1B2, hitsA2B1, hitsA2B2, this->searchType);
+ }
+ }
+ }
+ // "A" is being passed via STDIN.
+ else {
+
+ BEDPE a;
+ // process each entry in A
+ while (getline(cin, bedLine)) {
+
+ // split the current line into ditinct fields
+ vector<string> bedFields;
+ Tokenize(bedLine,bedFields);
+
+ lineNum++;
+
+ // find the overlaps with B if it's a valid BED entry.
+ if (bedA->parseBedPELine(a, bedFields, lineNum)) {
+ vector<BED> hitsA1B1, hitsA1B2, hitsA2B1, hitsA2B2;
+
+ FindOverlaps(a, hitsA1B1, hitsA1B2, hitsA2B1, hitsA2B2, this->searchType);
+ }
+ }
+ }
+}
+// END IntersectPE
+
+
+
+void PairToPair::FindOverlaps(BEDPE &a, vector<BED> &hitsA1B1, vector<BED> &hitsA1B2,
+ vector<BED> &hitsA2B1, vector<BED> &hitsA2B2, string &type) {
+
+ // list of hits on each end of BEDPE
+ // that exceed the requested overlap fraction
+ vector<BED> qualityHitsA1B1;
+ vector<BED> qualityHitsA1B2;
+ vector<BED> qualityHitsA2B1;
+ vector<BED> qualityHitsA2B2;
+
+ // count of hits on each end of BEDPE
+ // that exceed the requested overlap fraction
+ int numOverlapsA1B1 = 0;
+ int numOverlapsA1B2 = 0;
+ int numOverlapsA2B1 = 0;
+ int numOverlapsA2B2 = 0;
+
+
+ // Find the _potential_ hits between each end of A and B
+ bedB->binKeeperFind(bedB->bedMapEnd1[a.chrom1], a.start1, a.end1, hitsA1B1); // hits between A1 to B1
+ bedB->binKeeperFind(bedB->bedMapEnd1[a.chrom2], a.start2, a.end2, hitsA2B1); // hits between A2 to B1
+ bedB->binKeeperFind(bedB->bedMapEnd2[a.chrom1], a.start1, a.end1, hitsA1B2); // hits between A1 to B2
+ bedB->binKeeperFind(bedB->bedMapEnd2[a.chrom2], a.start2, a.end2, hitsA2B2); // hits between A2 to B2
+
+
+ // Now, reduce to the set of hits on each end of A and B
+ // that meet the required overlap fraction.
+ FindQualityHitsBetweenEnds(a, 1, hitsA1B1, qualityHitsA1B1, numOverlapsA1B1);
+ FindQualityHitsBetweenEnds(a, 1, hitsA1B2, qualityHitsA1B2, numOverlapsA1B2);
+ FindQualityHitsBetweenEnds(a, 2, hitsA2B1, qualityHitsA2B1, numOverlapsA2B1);
+ FindQualityHitsBetweenEnds(a, 2, hitsA2B2, qualityHitsA2B2, numOverlapsA2B2);
+
+
+ int matchCount1 = 0;
+ int matchCount2 = 0;
+ if ((numOverlapsA1B1 > 0) && (numOverlapsA2B2 > 0)) {
+
+ map<unsigned int, vector<BED>, less<unsigned int> > hitsMap;
+
+ for (vector<BED>::iterator h = qualityHitsA1B1.begin(); h != qualityHitsA1B1.end(); ++h) {
+ hitsMap[h->count].push_back(*h);
+ }
+ for (vector<BED>::iterator h = qualityHitsA2B2.begin(); h != qualityHitsA2B2.end(); ++h) {
+ hitsMap[h->count].push_back(*h);
+ }
+ for (map<unsigned int, vector<BED>, less<unsigned int> >::iterator m = hitsMap.begin(); m != hitsMap.end(); ++m) {
+
+ if (m->second.size() == 2) {
+ BED b1 = m->second[0];
+ BED b2 = m->second[1];
+
+ matchCount1++;
+
+ if (this->searchType != "neither") {
+ bedA->reportBedPETab(a);
+ printf("%s\t%d\t%d\t%s\t%d\t%d\t%s\t%s\t%s\t%s\n", b1.chrom.c_str(), b1.start, b1.end,
+ b2.chrom.c_str(), b2.start, b2.end,
+ b1.name.c_str(), b1.score.c_str(), b1.strand.c_str(), b2.strand.c_str());
+ }
+ }
+ }
+ }
+
+ if ((numOverlapsA1B2 > 0) && (numOverlapsA2B1 > 0)) {
+
+ map<unsigned int, vector<BED>, less<int> > hitsMap;
+
+ for (vector<BED>::iterator h = qualityHitsA2B1.begin(); h != qualityHitsA2B1.end(); ++h) {
+ hitsMap[h->count].push_back(*h);
+ }
+ for (vector<BED>::iterator h = qualityHitsA1B2.begin(); h != qualityHitsA1B2.end(); ++h) {
+ hitsMap[h->count].push_back(*h);
+ }
+ for (map<unsigned int, vector<BED>, less<unsigned int> >::iterator m = hitsMap.begin(); m != hitsMap.end(); ++m) {
+ if (m->second.size() == 2) {
+
+ BED b1 = m->second[0];
+ BED b2 = m->second[1];
+
+ matchCount2++;
+
+ if (this->searchType != "neither") {
+ bedA->reportBedPETab(a);
+ printf("%s\t%d\t%d\t%s\t%d\t%d\t%s\t%s\t%s\t%s\n", b1.chrom.c_str(), b1.start, b1.end,
+ b2.chrom.c_str(), b2.start, b2.end,
+ b1.name.c_str(), b1.score.c_str(), b1.strand.c_str(), b2.strand.c_str());
+ }
+ }
+ }
+ }
+
+ if ((matchCount1 == 0) && (matchCount2 == 0) && (this->searchType == "neither")) {
+ bedA->reportBedPENewLine(a);
+ }
+}
+
+
+
+void PairToPair::FindQualityHitsBetweenEnds(BEDPE a, int end, vector<BED> &hits, vector<BED> &qualityHits, int &numOverlaps) {
+
+ if (end == 1) {
+ for (vector<BED>::iterator h = hits.begin(); h != hits.end(); ++h) {
+ int s = max(a.start1, h->start);
+ int e = min(a.end1, h->end);
+ if (s < e) {
+ // is there enough overlap (default ~ 1bp)
+ if ( ((float)(e-s) / (float)(a.end1 - a.start1)) >= this->overlapFraction ) {
+ numOverlaps++;
+ qualityHits.push_back(*h);
+ }
+ }
+ }
+ }
+ else if (end == 2) {
+ for (vector<BED>::iterator h = hits.begin(); h != hits.end(); ++h) {
+ int s = max(a.start2, h->start);
+ int e = min(a.end2, h->end);
+ if (s < e) {
+ // is there enough overlap (default ~ 1bp)
+ if ( ((float)(e-s) / (float)(a.end2 - a.start2)) >= this->overlapFraction ) {
+ numOverlaps++;
+ qualityHits.push_back(*h);
+ }
+ }
+ }
+ }
+}
diff --git a/src/pairToPair/pairToPair.h b/src/pairToPair/pairToPair.h
new file mode 100755
index 0000000000000000000000000000000000000000..da3d25e4f96332a256cc173ce225c6c6af58200d
--- /dev/null
+++ b/src/pairToPair/pairToPair.h
@@ -0,0 +1,46 @@
+#ifndef PAIRTOPAIR_H
+#define PAIRTOPAIR_H
+
+#include "bedFile.h"
+#include "bedFilePE.h"
+#include <vector>
+#include <iostream>
+#include <fstream>
+
+using namespace std;
+
+//************************************************
+// Class methods and elements
+//************************************************
+class PairToPair {
+
+public:
+
+ // constructor
+ PairToPair(string &, string &, float &, string &);
+
+ // destructor
+ ~PairToPair(void);
+
+ void IntersectPairs ();
+
+ void FindOverlaps(BEDPE &, vector<BED> &, vector<BED> &, vector<BED> &, vector<BED> &, string &);
+
+ void FindQualityHitsBetweenEnds(BEDPE, int, vector<BED> &, vector<BED> &, int &);
+
+private:
+
+ string bedAFilePE;
+ string bedBFilePE;
+
+ float overlapFraction;
+ string searchType;
+
+ // instance of a paired-end bed file class.
+ BedFilePE *bedA;
+
+ // instance of a bed file class.
+ BedFilePE *bedB;
+};
+
+#endif /* PAIRTOPAIR_H */
diff --git a/src/pairToPair/pairToPairMain.cpp b/src/pairToPair/pairToPairMain.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..50d029a30f55633f8e58e80d71ba7ae167f7b16f
--- /dev/null
+++ b/src/pairToPair/pairToPairMain.cpp
@@ -0,0 +1,148 @@
+#include "pairToPair.h"
+#include "version.h"
+
+using namespace std;
+
+// define our program name
+#define PROGRAM_NAME "peIntersectBed"
+
+// 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;
+ string searchType = "both";
+
+ // flags to track parameters
+ bool haveBedA = false;
+ bool haveBedB = false;
+ bool haveSearchType = false;
+ bool haveFraction = false;
+
+ // 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;
+ 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("-type", 5, parameterLength)) {
+ if ((i+1) < argc) {
+ haveSearchType = true;
+ searchType = argv[i + 1];
+ }
+ i++;
+ }
+ else if(PARAMETER_CHECK("-f", 2, parameterLength)) {
+ haveFraction = true;
+ overlapFraction = atof(argv[i + 1]);
+ i++;
+ }
+ 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 (haveSearchType && (searchType != "neither") && (searchType != "both")) {
+ cerr << endl << "*****" << endl << "*****ERROR: Request \"both\" or \"neither\"" << endl << "*****" << endl;
+ showHelp = true;
+ }
+
+ if (!showHelp) {
+
+ PairToPair *bi = new PairToPair(bedAFile, bedBFile, overlapFraction, searchType);
+ bi->IntersectPairs();
+ return 0;
+ }
+ else {
+ ShowHelp();
+ }
+}
+
+void ShowHelp(void) {
+
+ cerr << "===============================================" << endl;
+ cerr << " " <<PROGRAM_NAME << " v" << VERSION << endl ;
+ cerr << " Aaron Quinlan, Ph.D. (aaronquinlan@gmail.com) " << endl ;
+ cerr << " Hall Laboratory, University of Virginia" << endl;
+ cerr << "===============================================" << endl << endl;
+ cerr << "Description: Report overlaps between a.bedpe and b.bedpe." << endl << endl;
+
+ cerr << "Usage: " << PROGRAM_NAME << " [OPTIONS] -a <a.bedpe> -b <b.bedpe>" << endl << endl;
+
+ cerr << "OPTIONS: " << endl;
+ cerr << "\t" << "-f\t\t" << "Minimum overlap req'd as fraction of a.bed (e.g. 0.05)." << endl << "\t\t\tDefault is 1E-9 (effectively 1bp)." << endl << endl;
+ cerr << "\t" << "-type \t\t" << "neither\t\t\tReport overlaps if _neither_ end of BEDPE (A) overlaps BEDPE B." << endl;
+ cerr << "\t\t\tboth\t\t\tReport overlaps if _both_ ends of BEDPE (A) overlap BEDPE B." << endl;
+
+
+ cerr << "NOTES: " << endl;
+ cerr << "\t" << "-i stdin\t" << "Allows BEDPE file A to be read from stdin. E.g.: cat a.bedpe | peIntersectBed -a stdin -b B.bed" << endl << endl;
+ cerr << "\t***Only 10 columns BEDPE formats allowed. See below).***"<< endl << endl;
+
+ cerr << "BEDPE FORMAT (tab-delimited): " << endl;
+ cerr << "\t" << "1. chrom for end 1" << endl;
+ cerr << "\t" << "2. start for end 1" << endl;
+ cerr << "\t" << "3. end for end 1" << endl;
+ cerr << "\t" << "4. chrom for end 2" << endl;
+ cerr << "\t" << "5. start for end 2" << endl;
+ cerr << "\t" << "6. end for end 2" << endl;
+ cerr << "\t" << "7. name" << endl;
+ cerr << "\t" << "8. score" << endl;
+ cerr << "\t" << "9. strand for end 1" << endl;
+ cerr << "\t" << "10. strand for end 2" << endl;
+ cerr << "\t" << "Note: Strands for each end must be provided if you choose to include strand information." << endl << endl;
+
+ cerr << "Example BEDPE record:" << endl;
+ cerr << "chr1\t1000000\t1000100\tchr1\t1001000\t1001100\tsomename\t100\t+\t-" << endl;
+
+ // end the program here
+ exit(1);
+
+}
diff --git a/src/peIntersectBed/peIntersectMain.cpp b/src/peIntersectBed/peIntersectMain.cpp
index 0dc22802b66872f0ad303ebbe8a2b60ca339d9f1..8931fa899eaf0e37543e04c34d0136d52cd2963e 100755
--- a/src/peIntersectBed/peIntersectMain.cpp
+++ b/src/peIntersectBed/peIntersectMain.cpp
@@ -113,9 +113,9 @@ void ShowHelp(void) {
cerr << " Aaron Quinlan, Ph.D. (aaronquinlan@gmail.com) " << endl ;
cerr << " Hall Laboratory, University of Virginia" << endl;
cerr << "===============================================" << endl << endl;
- cerr << "Description: Report overlaps between a.bed and b.bed." << endl << endl;
+ cerr << "Description: Report overlaps between a.bedpe and b.bed." << endl << endl;
- cerr << "Usage: " << PROGRAM_NAME << " [OPTIONS] -a <a.bed> -b <b.bed>" << endl << endl;
+ cerr << "Usage: " << PROGRAM_NAME << " [OPTIONS] -a <a.bedpe> -b <b.bed>" << endl << endl;
cerr << "OPTIONS: " << endl;
//cerr << "\t" << "-u\t\t" << "Write ORIGINAL a.bed entry ONCE if ANY overlap with B.bed." << endl << "\t\t\tIn other words, just report the fact >=1 hit was found." << endl << endl;
diff --git a/src/utils/bedFilePE/Makefile b/src/utils/bedFilePE/Makefile
index b99a1600cc3bc20ea33f65195e5cda83346e4e8f..f3853d54ce6f93d53de24711103e276a8e70dae8 100755
--- a/src/utils/bedFilePE/Makefile
+++ b/src/utils/bedFilePE/Makefile
@@ -7,7 +7,7 @@ UTILITIES_DIR = ../../utils/
# -------------------
# define our includes
# -------------------
-INCLUDES = -I$(UTILITIES_DIR)/lineFileUtilities/
+INCLUDES = -I$(UTILITIES_DIR)/lineFileUtilities/ -I$(UTILITIES_DIR)/bedFile/
# ----------------------------------
# define our source and object files
diff --git a/src/utils/bedFilePE/bedFilePE.cpp b/src/utils/bedFilePE/bedFilePE.cpp
index 937150dc9350060272c8589e65d11aafecebb87b..2cbbe798f3682c8a118e35390719d99cba361c2f 100755
--- a/src/utils/bedFilePE/bedFilePE.cpp
+++ b/src/utils/bedFilePE/bedFilePE.cpp
@@ -14,6 +14,74 @@
#include "bedFilePE.h"
+static int binOffsetsExtended[] =
+ {4096+512+64+8+1, 512+64+8+1, 64+8+1, 8+1, 1, 0};
+
+
+#define _binFirstShift 17 /* How much to shift to get to finest bin. */
+#define _binNextShift 3 /* How much to shift to get to next larger bin. */
+
+
+/*
+ NOTE: Taken ~verbatim from kent source.
+ Return a list of all items in binKeeper that intersect range.
+
+ Free this list with slFreeList.
+*/
+void BedFilePE::binKeeperFind(map<int, vector<BED>, std::less<int> > &bk, const int start, const int end, vector<BED> &hits) {
+
+ int startBin, endBin;
+ startBin = (start >>_binFirstShift);
+ endBin = ((end-1) >>_binFirstShift);
+
+ for (int i = 0; i < 6; ++i) {
+ int offset = binOffsetsExtended[i];
+
+ for (int j = (startBin+offset); j <= (endBin+offset); ++j) {
+ for (vector<BED>::const_iterator el = bk[j].begin(); el != bk[j].end(); ++el) {
+
+ if (overlaps(el->start, el->end, start, end) > 0) {
+ hits.push_back(*el);
+ }
+
+ }
+ }
+ startBin >>= _binNextShift;
+ endBin >>= _binNextShift;
+ }
+}
+
+
+
+static int getBin(int start, int end)
+/*
+ NOTE: Taken ~verbatim from kent source.
+
+ Given start,end in chromosome coordinates assign it
+ * a bin. There's a bin for each 128k segment, for each
+ * 1M segment, for each 8M segment, for each 64M segment,
+ * and for each chromosome (which is assumed to be less than
+ * 512M.) A range goes into the smallest bin it will fit in. */
+{
+ int startBin = start;
+ int endBin = end-1;
+ startBin >>= _binFirstShift;
+ endBin >>= _binFirstShift;
+
+ for (int i=0; i<6; ++i) {
+ if (startBin == endBin) {
+ return binOffsetsExtended[i] + startBin;
+ }
+ startBin >>= _binNextShift;
+ endBin >>= _binNextShift;
+ }
+
+ cerr << "start " << start << ", end " << end << " out of range in findBin (max is 512M)" << endl;
+ return 0;
+}
+
+
+
//***********************************************
// Common functions
//***********************************************
@@ -46,14 +114,14 @@ void BedFilePE::reportBedPETab(const BEDPE &a) {
a.name.c_str());
}
else if (this->bedType == 8) {
- printf("%s\t%d\t%d\t%s\t%d\t%d\t%s\t%d\t", a.chrom1.c_str(), a.start1, a.end1,
+ printf("%s\t%d\t%d\t%s\t%d\t%d\t%s\t%s\t", a.chrom1.c_str(), a.start1, a.end1,
a.chrom2.c_str(), a.start2, a.end2,
- a.name.c_str(), a.score);
+ a.name.c_str(), a.score.c_str());
}
else if (this->bedType == 10) {
- printf("%s\t%d\t%d\t%s\t%d\t%d\t%s\t%d\t%s\t%s\t", a.chrom1.c_str(), a.start1, a.end1,
+ printf("%s\t%d\t%d\t%s\t%d\t%d\t%s\t%s\t%s\t%s\t", a.chrom1.c_str(), a.start1, a.end1,
a.chrom2.c_str(), a.start2, a.end2,
- a.name.c_str(), a.score, a.strand1.c_str(), a.strand2.c_str());
+ a.name.c_str(), a.score.c_str(), a.strand1.c_str(), a.strand2.c_str());
}
}
@@ -77,14 +145,14 @@ void BedFilePE::reportBedPENewLine(const BEDPE &a) {
a.name.c_str());
}
else if (this->bedType == 8) {
- printf("%s\t%d\t%d\t%s\t%d\t%d\t%s\t%d\n", a.chrom1.c_str(), a.start1, a.end1,
+ printf("%s\t%d\t%d\t%s\t%d\t%d\t%s\t%s\n", a.chrom1.c_str(), a.start1, a.end1,
a.chrom2.c_str(), a.start2, a.end2,
- a.name.c_str(), a.score);
+ a.name.c_str(), a.score.c_str());
}
else if (this->bedType == 10) {
- printf("%s\t%d\t%d\t%s\t%d\t%d\t%s\t%d\t%s\t%s\n", a.chrom1.c_str(), a.start1, a.end1,
+ printf("%s\t%d\t%d\t%s\t%d\t%d\t%s\t%s\t%s\t%s\n", a.chrom1.c_str(), a.start1, a.end1,
a.chrom2.c_str(), a.start2, a.end2,
- a.name.c_str(), a.score, a.strand1.c_str(), a.strand2.c_str());
+ a.name.c_str(), a.score.c_str(), a.strand1.c_str(), a.strand2.c_str());
}
}
@@ -256,3 +324,107 @@ bool BedFilePE::parseBedPELine (BEDPE &bed, const vector<string> &lineVector, co
}
+void BedFilePE::loadBedPEFileIntoMap() {
+
+ // open the BEDPE file for reading
+ ifstream bed(bedFile.c_str(), ios::in);
+ if ( !bed ) {
+ cerr << "Error: The requested bed file (" <<bedFile << ") could not be opened. Exiting!" << endl;
+ exit (1);
+ }
+
+ string bedLine;
+ BEDPE bedpeEntry;
+ BED bedEntry1, bedEntry2;
+
+ unsigned int lineNum = 0;
+ int bin1, bin2;
+ vector<string> bedFields; // vector of strings for each column in BED file.
+ bedFields.reserve(10); // reserve space for worst case (BED 6)
+
+ while (getline(bed, bedLine)) {
+ lineNum++;
+
+ if (lineNum > 1) {
+
+ Tokenize(bedLine,bedFields);
+
+ if (parseBedPELine(bedpeEntry, bedFields, lineNum)) {
+
+ // separate the BEDPE entry into separate
+ // BED entries
+ splitBedPEIntoBeds(bedpeEntry, lineNum, bedEntry1, bedEntry2);
+
+ // load end1 into a UCSC bin map
+ bin1 = getBin(bedEntry1.start, bedEntry1.end);
+ this->bedMapEnd1[bedEntry1.chrom][bin1].push_back(bedEntry1);
+
+ // load end2 into a UCSC bin map
+ bin2 = getBin(bedEntry2.start, bedEntry2.end);
+ this->bedMapEnd2[bedEntry2.chrom][bin2].push_back(bedEntry2);
+ }
+ bedFields.clear();
+
+ }
+ else {
+ if ((bedLine.find("track") != string::npos) || (bedLine.find("browser") != string::npos)) {
+ continue;
+ }
+ else {
+ Tokenize(bedLine,bedFields);
+
+ if (parseBedPELine(bedpeEntry, bedFields, lineNum)) {
+ // separate the BEDPE entry into separate
+ // BED entries
+ splitBedPEIntoBeds(bedpeEntry, lineNum, bedEntry1, bedEntry2);
+
+ // load end1 into a UCSC bin map
+ bin1 = getBin(bedEntry1.start, bedEntry1.end);
+ this->bedMapEnd1[bedEntry1.chrom][bin1].push_back(bedEntry1);
+
+ // load end2 into a UCSC bin map
+ bin2 = getBin(bedEntry2.start, bedEntry2.end);
+ this->bedMapEnd2[bedEntry2.chrom][bin2].push_back(bedEntry2);
+ }
+ bedFields.clear();
+ }
+ }
+ }
+}
+
+
+void BedFilePE::splitBedPEIntoBeds(const BEDPE &bedpeEntry, unsigned int lineNum, BED &bedEntry1, BED &bedEntry2) {
+
+ /*
+ Split the BEDPE entry into separate BED entries
+
+ NOTE: I am using a trick here where I store
+ the lineNum of the BEDPE from the original file
+ in the "count" column. This allows me to later
+ resolve whether the hits found on both ends of BEDPE A
+ came from the same entry in BEDPE B. Tracking by "name"
+ alone with fail when there are multiple mappings for a given
+ read-pair.
+ */
+
+ bedEntry1.chrom = bedpeEntry.chrom1;
+ bedEntry1.start = bedpeEntry.start1;
+ bedEntry1.end = bedpeEntry.end1;
+ bedEntry1.name = bedpeEntry.name;
+ bedEntry1.score = bedpeEntry.score;
+ bedEntry1.strand = bedpeEntry.strand1;
+ bedEntry1.count = lineNum;
+ bedEntry1.minOverlapStart = INT_MAX;
+
+ bedEntry2.chrom = bedpeEntry.chrom2;
+ bedEntry2.start = bedpeEntry.start2;
+ bedEntry2.end = bedpeEntry.end2;
+ bedEntry2.name = bedpeEntry.name;
+ bedEntry2.score = bedpeEntry.score;
+ bedEntry2.strand = bedpeEntry.strand2;
+ bedEntry2.count = lineNum;
+ bedEntry2.minOverlapStart = INT_MAX;
+}
+
+
+
diff --git a/src/utils/bedFilePE/bedFilePE.h b/src/utils/bedFilePE/bedFilePE.h
index 9d26020d926f20627d4332a52ff88c5979f4033b..a48ce77d8b16e6d96f1699614f87565180b1417a 100755
--- a/src/utils/bedFilePE/bedFilePE.h
+++ b/src/utils/bedFilePE/bedFilePE.h
@@ -9,6 +9,8 @@
#include <sstream>
#include <cstring>
#include <algorithm>
+#include "bedFile.h"
+#include "lineFileUtilities.h"
using namespace std;
@@ -28,7 +30,7 @@ struct BEDPE {
int end2;
string name;
- int score;
+ string score;
string strand1;
string strand2;
@@ -36,6 +38,7 @@ struct BEDPE {
+
//************************************************
// BedFile Class methods and elements
//************************************************
@@ -53,11 +56,19 @@ public:
bool parseBedPELine (BEDPE &, const vector<string> &, const int &);
void reportBedPETab(const BEDPE &);
void reportBedPENewLine(const BEDPE &);
-
+ void loadBedPEFileIntoMap();
+ void splitBedPEIntoBeds(const BEDPE &, unsigned int, BED &, BED &);
+ void binKeeperFind(map<int, vector<BED>,
+ std::less<int> > &, const int,
+ const int, vector<BED> &);
+
string bedFile;
unsigned int bedType;
-
+
+ masterBedMap bedMapEnd1;
+ masterBedMap bedMapEnd2;
+
private:
// none
};