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

Issue #50. SlopBed now adds percent-based slop to features with the -pct option.

parent 0603eb52
No related branches found
No related tags found
No related merge requests found
...@@ -7,24 +7,25 @@ ...@@ -7,24 +7,25 @@
University of Virginia University of Virginia
aaronquinlan@gmail.com aaronquinlan@gmail.com
Licenced under the GNU General Public License 2.0 license. Licensed under the GNU General Public License 2.0 license.
******************************************************************************/ ******************************************************************************/
#include "lineFileUtilities.h" #include "lineFileUtilities.h"
#include "slopBed.h" #include "slopBed.h"
BedSlop::BedSlop(string &bedFile, string &genomeFile, bool &forceStrand, int &leftSlop, int &rightSlop) { BedSlop::BedSlop(string &bedFile, string &genomeFile, bool forceStrand, float leftSlop, float rightSlop, bool fractional) {
_bedFile = bedFile; _bedFile = bedFile;
_genomeFile = genomeFile; _genomeFile = genomeFile;
_forceStrand = forceStrand; _forceStrand = forceStrand;
_leftSlop = leftSlop;
_leftSlop = leftSlop; _rightSlop = rightSlop;
_rightSlop = rightSlop; _fractional = fractional;
_bed = new BedFile(bedFile); _bed = new BedFile(bedFile);
_genome = new GenomeFile(genomeFile); _genome = new GenomeFile(genomeFile);
// get going, slop it up.
SlopBed(); SlopBed();
} }
...@@ -44,7 +45,14 @@ void BedSlop::SlopBed() { ...@@ -44,7 +45,14 @@ void BedSlop::SlopBed() {
bedStatus = _bed->GetNextBed(bedEntry, lineNum); bedStatus = _bed->GetNextBed(bedEntry, lineNum);
while (bedStatus != BED_INVALID) { while (bedStatus != BED_INVALID) {
if (bedStatus == BED_VALID) { if (bedStatus == BED_VALID) {
AddSlop(bedEntry); if (_fractional == false) {
AddSlop(bedEntry, _leftSlop, _rightSlop);
}
else {
int leftSlop = (int) (_leftSlop * bedEntry.size());
int rightSlop = (int) (_rightSlop * bedEntry.size());
AddSlop(bedEntry, leftSlop, rightSlop);
}
_bed->reportBedNewLine(bedEntry); _bed->reportBedNewLine(bedEntry);
bedEntry = nullBed; bedEntry = nullBed;
} }
...@@ -54,7 +62,7 @@ void BedSlop::SlopBed() { ...@@ -54,7 +62,7 @@ void BedSlop::SlopBed() {
} }
void BedSlop::AddSlop(BED &bed) { void BedSlop::AddSlop(BED &bed, int leftSlop, int rightSlop) {
// special handling if the BED entry is on the negative // special handling if the BED entry is on the negative
// strand and the user cares about strandedness. // strand and the user cares about strandedness.
...@@ -62,20 +70,20 @@ void BedSlop::AddSlop(BED &bed) { ...@@ -62,20 +70,20 @@ void BedSlop::AddSlop(BED &bed) {
if ( (_forceStrand) && (bed.strand == "-") ) { if ( (_forceStrand) && (bed.strand == "-") ) {
// inspect the start // inspect the start
if ( (static_cast<int>(bed.start) - _rightSlop) > 0 ) bed.start -= _rightSlop; if ( (static_cast<int>(bed.start) - rightSlop) > 0 ) bed.start -= rightSlop;
else bed.start = 0; else bed.start = 0;
// inspect the start // inspect the start
if ( (static_cast<int>(bed.end) + _leftSlop) <= static_cast<int>(chromSize)) bed.end += _leftSlop; if ( (static_cast<int>(bed.end) + leftSlop) <= static_cast<int>(chromSize)) bed.end += leftSlop;
else bed.end = chromSize; else bed.end = chromSize;
} }
else { else {
// inspect the start // inspect the start
if ( (static_cast<int>(bed.start) - _leftSlop) > 0) bed.start -= _leftSlop; if ( (static_cast<int>(bed.start) - leftSlop) > 0) bed.start -= leftSlop;
else bed.start = 0; else bed.start = 0;
// inspect the end // inspect the end
if ( (static_cast<int>(bed.end) + _rightSlop) <= static_cast<int>(chromSize)) bed.end += _rightSlop; if ( (static_cast<int>(bed.end) + rightSlop) <= static_cast<int>(chromSize)) bed.end += rightSlop;
else bed.end = chromSize; else bed.end = chromSize;
} }
} }
......
...@@ -30,7 +30,7 @@ class BedSlop { ...@@ -30,7 +30,7 @@ class BedSlop {
public: public:
// constructor // constructor
BedSlop(string &bedFile, string &genomeFile, bool &forceStrand, int &leftSlop, int &rightSlop) ; BedSlop(string &bedFile, string &genomeFile, bool forceStrand, float leftSlop, float rightSlop, bool fractional);
// destructor // destructor
~BedSlop(void); ~BedSlop(void);
...@@ -42,9 +42,10 @@ private: ...@@ -42,9 +42,10 @@ private:
string _bedFile; string _bedFile;
string _genomeFile; string _genomeFile;
bool _forceStrand; bool _forceStrand;
int _leftSlop; float _leftSlop;
int _rightSlop; float _rightSlop;
bool _fractional;
BedFile *_bed; BedFile *_bed;
GenomeFile *_genome; GenomeFile *_genome;
...@@ -54,5 +55,5 @@ private: ...@@ -54,5 +55,5 @@ private:
void SlopBed(); void SlopBed();
// method to add requested "slop" to a single BED entry // method to add requested "slop" to a single BED entry
void AddSlop(BED &bed); void AddSlop(BED &bed, int leftSlop, int rightSlop);
}; };
...@@ -40,8 +40,9 @@ int main(int argc, char* argv[]) { ...@@ -40,8 +40,9 @@ int main(int argc, char* argv[]) {
bool haveBoth = false; bool haveBoth = false;
bool forceStrand = false; bool forceStrand = false;
int leftSlop = 0; float leftSlop = 0.0;
int rightSlop = 0; float rightSlop = 0.0;
bool fractional = false;
for(int i = 1; i < argc; i++) { for(int i = 1; i < argc; i++) {
int parameterLength = (int)strlen(argv[i]); int parameterLength = (int)strlen(argv[i]);
...@@ -75,28 +76,31 @@ int main(int argc, char* argv[]) { ...@@ -75,28 +76,31 @@ int main(int argc, char* argv[]) {
else if(PARAMETER_CHECK("-l", 2, parameterLength)) { else if(PARAMETER_CHECK("-l", 2, parameterLength)) {
if ((i+1) < argc) { if ((i+1) < argc) {
haveLeft = true; haveLeft = true;
leftSlop = atoi(argv[i + 1]); leftSlop = atof(argv[i + 1]);
i++; i++;
} }
} }
else if(PARAMETER_CHECK("-r", 2, parameterLength)) { else if(PARAMETER_CHECK("-r", 2, parameterLength)) {
if ((i+1) < argc) { if ((i+1) < argc) {
haveRight = true; haveRight = true;
rightSlop = atoi(argv[i + 1]); rightSlop = atof(argv[i + 1]);
i++; i++;
} }
} }
else if(PARAMETER_CHECK("-b", 2, parameterLength)) { else if(PARAMETER_CHECK("-b", 2, parameterLength)) {
if ((i+1) < argc) { if ((i+1) < argc) {
haveBoth = true; haveBoth = true;
leftSlop = atoi(argv[i + 1]); leftSlop = atof(argv[i + 1]);
rightSlop = atoi(argv[i + 1]); rightSlop = atof(argv[i + 1]);
i++; i++;
} }
} }
else if(PARAMETER_CHECK("-s", 2, parameterLength)) { else if(PARAMETER_CHECK("-s", 2, parameterLength)) {
forceStrand = true; forceStrand = true;
} }
else if(PARAMETER_CHECK("-pct", 4, parameterLength)) {
fractional = true;
}
else { else {
cerr << endl << "*****ERROR: Unrecognized parameter: " << argv[i] << " *****" << endl << endl; cerr << endl << "*****ERROR: Unrecognized parameter: " << argv[i] << " *****" << endl << endl;
showHelp = true; showHelp = true;
...@@ -122,7 +126,7 @@ int main(int argc, char* argv[]) { ...@@ -122,7 +126,7 @@ int main(int argc, char* argv[]) {
} }
if (!showHelp) { if (!showHelp) {
BedSlop *bc = new BedSlop(bedFile, genomeFile, forceStrand, leftSlop, rightSlop); BedSlop *bc = new BedSlop(bedFile, genomeFile, forceStrand, leftSlop, rightSlop, fractional);
delete bc; delete bc;
return 0; return 0;
...@@ -144,18 +148,22 @@ void ShowHelp(void) { ...@@ -144,18 +148,22 @@ void ShowHelp(void) {
cerr << "Options: " << endl; cerr << "Options: " << endl;
cerr << "\t-b\t" << "Increase the BED/GFF/VCF entry by -b base pairs in each direction." << endl; cerr << "\t-b\t" << "Increase the BED/GFF/VCF entry by -b base pairs in each direction." << endl;
cerr << "\t\t- (Integer)" << endl; cerr << "\t\t- (Integer) or (Float, e.g. 0.1) if used with -pct." << endl << endl;
cerr << "\t-l\t" << "The number of base pairs to subtract from the start coordinate." << endl; cerr << "\t-l\t" << "The number of base pairs to subtract from the start coordinate." << endl;
cerr << "\t\t- (Integer)" << endl; cerr << "\t\t- (Integer) or (Float, e.g. 0.1) if used with -pct." << endl << endl;
cerr << "\t-r\t" << "The number of base pairs to add to the end coordinate." << endl; cerr << "\t-r\t" << "The number of base pairs to add to the end coordinate." << endl;
cerr << "\t\t- (Integer)" << endl; cerr << "\t\t- (Integer) or (Float, e.g. 0.1) if used with -pct." << endl << endl;
cerr << "\t-s\t" << "Define -l and -r based on strand." << endl; cerr << "\t-s\t" << "Define -l and -r based on strand." << endl;
cerr << "\t\tE.g. if used, -l 500 for a negative-stranded feature, " << endl; cerr << "\t\tE.g. if used, -l 500 for a negative-stranded feature, " << endl;
cerr << "\t\tit will add 500 bp downstream. Default = false." << endl << endl; cerr << "\t\tit will add 500 bp downstream. Default = false." << endl << endl;
cerr << "\t-pct\t" << "Define -l and -r as a fraction of the feature's length." << endl;
cerr << "\t\tE.g. if used on a 1000bp feature, -l 0.50, " << endl;
cerr << "\t\twill add 500 bp \"upstream\". Default = false." << endl << endl;
cerr << "Notes: " << endl; cerr << "Notes: " << endl;
cerr << "\t(1) Starts will be set to 0 if options would force it below 0." << endl; cerr << "\t(1) Starts will be set to 0 if options would force it below 0." << endl;
cerr << "\t(2) Ends will be set to the chromosome length if requested slop would" << endl; cerr << "\t(2) Ends will be set to the chromosome length if requested slop would" << endl;
......
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