Skip to content
Snippets Groups Projects
Commit 9713b9fb authored by Jake Biesinger's avatar Jake Biesinger
Browse files

User can specify one of three orientations (a,b,ref)

parent a65203c8
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/closestBed/closestBed.cpp
+ 31
− 53
View file @ 9713b9fb
......@@ -22,7 +22,8 @@ const int SLOPGROWTH = 2048000;
Constructor
*/
BedClosest::BedClosest(string &bedAFile, string &bedBFile, bool sameStrand, bool diffStrand,
string &tieMode, bool reportDistance, bool signDistance, bool ignoreOverlaps)
string &tieMode, bool reportDistance, bool signDistance, string &_strandedDistMode,
bool ignoreOverlaps)
: _bedAFile(bedAFile)
, _bedBFile(bedBFile)
, _tieMode(tieMode)
......@@ -30,6 +31,7 @@ BedClosest::BedClosest(string &bedAFile, string &bedBFile, bool sameStrand, bool
, _diffStrand(diffStrand)
, _reportDistance(reportDistance)
, _signDistance(signDistance)
, _strandedDistMode(_strandedDistMode)
, _ignoreOverlaps(ignoreOverlaps)
{
_bedA = new BedFile(_bedAFile);
......@@ -56,7 +58,8 @@ void BedClosest::FindWindowOverlaps(BED &a, vector<BED> &hits) {
CHRPOS aFudgeEnd;
int numOverlaps = 0;
vector<BED> closestB;
int32_t minDistance = INT_MAX;
CHRPOS minDistance = INT_MAX;
int32_t curDistance = INT_MAX;
vector<int32_t> distances;
// is there at least one feature in B on the same chrom
......@@ -103,73 +106,48 @@ void BedClosest::FindWindowOverlaps(BED &a, vector<BED> &hits) {
}
// the hit is to the "left" of A
else if (h->end <= a.start) {
if ((a.start - h->end) < abs(minDistance)) {
if (_signDistance) {
if (a.strand == "+") {
minDistance = h->end - a.start;
}
else {
minDistance = a.start - h->end;
}
curDistance = a.start - h->end;
if (_signDistance) {
if ((_strandedDistMode == "ref")
|| (_strandedDistMode == "a" && a.strand != "-")
|| (_strandedDistMode == "b" && h->strand == "-")) {
curDistance = -curDistance;
}
else {
minDistance = a.start - h->end;
}
}
if (abs(curDistance) < minDistance) {
minDistance = abs(curDistance);
closestB.clear();
closestB.push_back(*h);
distances.clear();
distances.push_back(minDistance);
distances.push_back(curDistance);
}
else if ((a.start - h->end) == abs(minDistance)) {
if (_signDistance) {
if (a.strand == "+") {
minDistance = h->end - a.start;
}
else {
minDistance = a.start - h->end;
}
}
else {
minDistance = h->end - a.start;
}
else if (abs(curDistance) == minDistance) {
minDistance = abs(curDistance);
closestB.push_back(*h);
distances.push_back(minDistance);
distances.push_back(curDistance);
}
}
// the hit is to the "right" of A
else if (h->start >= a.end) {
if ((h->start - a.end) < abs(minDistance)) {
if (_signDistance) {
if (a.strand == "+") {
minDistance = h->start - a.end;
}
else {
minDistance = a.end - h->start;
}
}
else {
minDistance = h->start - a.end;
curDistance = h->start - a.end;
if (_signDistance) {
if ((_strandedDistMode == "a" && a.strand == "-")
|| (_strandedDistMode == "b" && h->strand != "-")) {
curDistance = -curDistance;
}
}
if (abs(curDistance) < minDistance) {
closestB.clear();
closestB.push_back(*h);
distances.clear();
distances.push_back(minDistance);
distances.push_back(curDistance);
}
else if ((h->start - a.end) == abs(minDistance)) {
if (_signDistance) {
if (a.strand == "+") {
minDistance = h->start - a.end;
}
else {
minDistance = a.end - h->start;
}
}
else {
minDistance = h->start - a.end;
}
else if (abs(curDistance) == minDistance) {
closestB.push_back(*h);
distances.push_back(minDistance);
distances.push_back(curDistance);
}
}
}
......
src/closestBed/closestBed.h
+ 3
− 1
View file @ 9713b9fb
......@@ -29,7 +29,8 @@ public:
// constructor
BedClosest(string &bedAFile, string &bedBFile,
bool sameStrand, bool diffStrand, string &tieMode,
bool reportDistance, bool signDistance, bool ignoreOverlaps);
bool reportDistance, bool signDistance, string &strandedDistMode,
bool ignoreOverlaps);
// destructor
~BedClosest(void);
......@@ -47,6 +48,7 @@ private:
bool _diffStrand;
bool _reportDistance;
bool _signDistance;
string _strandedDistMode;
bool _ignoreOverlaps;
BedFile *_bedA, *_bedB;
......
src/closestBed/closestMain.cpp
+ 24
− 5
View file @ 9713b9fb
......@@ -32,6 +32,7 @@ int main(int argc, char* argv[]) {
string bedAFile;
string bedBFile;
string tieMode = "all";
string strandedDistMode = "ref";
bool haveBedA = false;
bool haveBedB = false;
......@@ -41,6 +42,7 @@ int main(int argc, char* argv[]) {
bool ignoreOverlaps = false;
bool reportDistance = false;
bool signDistance = false;
bool haveStrandedDistMode = false;
// check to see if we should print out some help
......@@ -86,8 +88,13 @@ int main(int argc, char* argv[]) {
reportDistance = true;
}
else if (PARAMETER_CHECK("-D", 2, parameterLength)) {
reportDistance = true;
signDistance = true;
if ((i+1) < argc) {
reportDistance = true;
signDistance = true;
haveStrandedDistMode = true;
strandedDistMode = argv[i + 1];
i++;
}
}
else if (PARAMETER_CHECK("-io", 3, parameterLength)) {
ignoreOverlaps = true;
......@@ -116,6 +123,12 @@ int main(int argc, char* argv[]) {
cerr << endl << "*****" << endl << "*****ERROR: Request \"all\" or \"first\" or \"last\" for Tie Mode (-t)" << endl << "*****" << endl;
showHelp = true;
}
if (haveStrandedDistMode && (strandedDistMode != "a") && (strandedDistMode != "b")
&& (strandedDistMode != "ref")) {
cerr << endl << "*****" << endl << "*****ERROR: Request \"a\" or \"b\" or \"ref\" for Stranded Distance Mode (-D)" << endl << "*****" << endl;
showHelp = true;
}
if (sameStrand && diffStrand) {
cerr << endl << "*****" << endl << "*****ERROR: Request either -s OR -S, not both." << endl << "*****" << endl;
......@@ -123,7 +136,7 @@ int main(int argc, char* argv[]) {
}
if (!showHelp) {
BedClosest *bc = new BedClosest(bedAFile, bedBFile, sameStrand, diffStrand, tieMode, reportDistance, signDistance, ignoreOverlaps);
BedClosest *bc = new BedClosest(bedAFile, bedBFile, sameStrand, diffStrand, tieMode, reportDistance, signDistance, strandedDistMode, ignoreOverlaps);
delete bc;
return 0;
}
......@@ -159,8 +172,14 @@ void ShowHelp(void) {
cerr << "\t-D\t" << "Like -d, report the closest feature in B, and its distance to A" << endl;
cerr << "\t\tas an extra column. Unlike -d, use negative distances to report" << endl;
cerr << "\t\tif feature B is upstream of A. Note: Upstream B features are" << endl;
cerr << "\"to the right\" of features on the - strand" << endl << endl;
cerr << "\t\tupstream features. You must specify which orientation defines \"upstream\"." << endl;
cerr << "\t\tThe options are:" << endl;
cerr << "\t\t- \"ref\" Report distance with respect to the reference genome. " << endl;
cerr << "\t\t B features with a lower (start, stop) are upstream" << endl;
cerr << "\t\t- \"a\" Report distance with respect to A." << endl;
cerr << "\t\t When A is on the - strand, \"upstream\" means B has a higher (start,stop)." << endl;
cerr << "\t\t- \"b\" Report distance with respect to B." << endl;
cerr << "\t\t When B is on the - strand, \"upstream\" means A has a higher (start,stop)." << endl << endl;
cerr << "\t-no\t" << "Ignore features in B that overlap A. That is, we want close, but " << endl;
cerr << "\t\tnot touching features only." << endl << endl;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment