From 832c4de2a9259aa0fb7f5e0df4d7eec5e40253e4 Mon Sep 17 00:00:00 2001
From: Aaron <aaronquinlan@gmail.com>
Date: Wed, 26 Oct 2011 21:40:51 -0400
Subject: [PATCH] nucBed now obeys strand'
---
src/nucBed/nucBed.cpp | 9 +++++++--
src/nucBed/nucBed.h | 5 ++++-
src/nucBed/nucBedMain.cpp | 7 ++++++-
3 files changed, 17 insertions(+), 4 deletions(-)
diff --git a/src/nucBed/nucBed.cpp b/src/nucBed/nucBed.cpp
index 49383112..5adaa2cf 100644
--- a/src/nucBed/nucBed.cpp
+++ b/src/nucBed/nucBed.cpp
@@ -13,13 +13,15 @@
#include "nucBed.h"
-NucBed::NucBed(string &dbFile, string &bedFile, bool printSeq, bool hasPattern, const string &pattern) {
+NucBed::NucBed(string &dbFile, string &bedFile, bool printSeq,
+ bool hasPattern, const string &pattern, bool forceStrand) {
_dbFile = dbFile;
_bedFile = bedFile;
_printSeq = printSeq;
_hasPattern = hasPattern;
_pattern = pattern;
+ _forceStrand = forceStrand;
_bed = new BedFile(_bedFile);
@@ -129,7 +131,10 @@ void NucBed::ProfileDNA() {
// grab the dna at this interval
int length = bed.end - bed.start;
// report the sequence's content
- ReportDnaProfile(bed, fr.getSubSequence(bed.chrom, bed.start, length), length);
+ string dna = fr.getSubSequence(bed.chrom, bed.start, length);
+ if ((_forceStrand == true) && (bed.strand == "-"))
+ reverseComplement(dna);
+ ReportDnaProfile(bed, dna, length);
bed = nullBed;
}
else
diff --git a/src/nucBed/nucBed.h b/src/nucBed/nucBed.h
index 323a8a49..2fe40ddf 100644
--- a/src/nucBed/nucBed.h
+++ b/src/nucBed/nucBed.h
@@ -29,7 +29,9 @@ class NucBed {
public:
// constructor
- NucBed(string &dbFile, string &bedFile, bool printSeq, bool hasPattern, const string &pattern);
+ NucBed(string &dbFile, string &bedFile, bool printSeq,
+ bool hasPattern, const string &pattern,
+ bool forceStrand);
// destructor
~NucBed(void);
@@ -42,6 +44,7 @@ private:
bool _printSeq;
bool _hasPattern;
string _pattern;
+ bool _forceStrand;
// instance of a bed file class.
BedFile *_bed;
diff --git a/src/nucBed/nucBedMain.cpp b/src/nucBed/nucBedMain.cpp
index 5fa45030..f92ddd88 100644
--- a/src/nucBed/nucBedMain.cpp
+++ b/src/nucBed/nucBedMain.cpp
@@ -39,6 +39,7 @@ int main(int argc, char* argv[]) {
bool haveBed = false;
bool printSeq = false;
bool hasPattern = false;
+ bool forceStrand = false;
// check to see if we should print out some help
if(argc <= 1) showHelp = true;
@@ -76,6 +77,9 @@ int main(int argc, char* argv[]) {
else if(PARAMETER_CHECK("-seq", 4, parameterLength)) {
printSeq = true;
}
+ else if(PARAMETER_CHECK("-s", 2, parameterLength)) {
+ forceStrand = true;
+ }
else if(PARAMETER_CHECK("-pattern", 8, parameterLength)) {
if ((i+1) < argc) {
hasPattern = true;
@@ -95,7 +99,7 @@ int main(int argc, char* argv[]) {
if (!showHelp) {
- NucBed *nuc = new NucBed(fastaDbFile, bedFile, printSeq, hasPattern, pattern);
+ NucBed *nuc = new NucBed(fastaDbFile, bedFile, printSeq, hasPattern, pattern, forceStrand);
delete nuc;
return 0;
@@ -118,6 +122,7 @@ void ShowHelp(void) {
cerr << "Options: " << endl;
cerr << "\t-fi\tInput FASTA file" << endl << endl;
cerr << "\t-bed\tBED/GFF/VCF file of ranges to extract from -fi" << endl << endl;
+ cerr << "\t-s\tProfile the sequence according to strand." << endl << endl;
cerr << "\t-seq\tPrint the extracted sequence" << endl << endl;
cerr << "\t-pattern\tReport the number of times a user-defined sequence is observed (case-insensitive)." << endl << endl;
--
GitLab