From 86f82fd1de975ac30dd54780a627905e8e5c7da1 Mon Sep 17 00:00:00 2001
From: Aaron <aaronquinlan@gmail.com>
Date: Wed, 28 Apr 2010 22:15:54 -0400
Subject: [PATCH] Added -bga mode to genomeCoverageBed. Thanks to Gordon
Assaf, CSHL.
---
src/genomeCoverageBed/genomeCoverageBed.cpp | 43 +++++++++-----------
src/genomeCoverageBed/genomeCoverageBed.h | 3 +-
src/genomeCoverageBed/genomeCoverageMain.cpp | 39 ++++++++++++------
3 files changed, 48 insertions(+), 37 deletions(-)
diff --git a/src/genomeCoverageBed/genomeCoverageBed.cpp b/src/genomeCoverageBed/genomeCoverageBed.cpp
index c0563a30..a506e183 100755
--- a/src/genomeCoverageBed/genomeCoverageBed.cpp
+++ b/src/genomeCoverageBed/genomeCoverageBed.cpp
@@ -14,15 +14,16 @@
BedGenomeCoverage::BedGenomeCoverage(string bedFile, string genomeFile, bool eachBase,
- bool startSites, bool bedGraph, int max, bool bamInput) {
-
- _bedFile = bedFile;
- _genomeFile = genomeFile;
- _eachBase = eachBase;
- _startSites = startSites;
- _bedGraph = bedGraph;
- _max = max;
- _bamInput = bamInput;
+ bool startSites, bool bedGraph, bool bedGraphAll, int max, bool bamInput) {
+
+ _bedFile = bedFile;
+ _genomeFile = genomeFile;
+ _eachBase = eachBase;
+ _startSites = startSites;
+ _bedGraph = bedGraph;
+ _bedGraphAll = bedGraphAll;
+ _max = max;
+ _bamInput = bamInput;
_bed = new BedFile(bedFile);
_genome = new GenomeFile(genomeFile);
@@ -124,10 +125,9 @@ void BedGenomeCoverage::CoverageBed() {
}
_bed->Close();
- cout << "yep\n";
// process the results of the last chromosome.
ReportChromCoverage(chromCov, currChromSize, currChrom, chromDepthHist);
- if (_eachBase == false && _bedGraph == false) {
+ if (_eachBase == false && _bedGraph == false && _bedGraphAll == false) {
ReportGenomeCoverage(chromDepthHist);
}
}
@@ -234,7 +234,7 @@ void BedGenomeCoverage::ReportChromCoverage(const vector<DEPTH> &chromCov, const
depth = depth - chromCov[pos].ends;
}
}
- else if (_bedGraph) {
+ else if (_bedGraph == true || _bedGraphAll == true) {
ReportChromCoverageBedGraph(chromCov, chromSize, chrom);
}
else {
@@ -323,16 +323,12 @@ void BedGenomeCoverage::ReportChromCoverageBedGraph(const vector<DEPTH> &chromCo
for (int pos = 0; pos < chromSize; pos++) {
depth += chromCov[pos].starts;
- if (depth == 0 && lastDepth != -1) {
- // We've found a new block of zero coverage, so report
- // the previous block of non-zero coverage.
- cout << chrom << "\t" << lastStart << "\t" << pos << "\t" << lastDepth << endl;
- lastDepth = -1;
- lastStart = -1;
- }
- else if (depth > 0 && depth != lastDepth) {
+ if (depth != lastDepth) {
// Coverage depth has changed, print the last interval coverage (if any)
- if (lastDepth != -1) {
+ // Print if:
+ // (1) depth>0 (the default running mode),
+ // (2) depth==0 and the user requested to print zero covered regions (_bedGraphAll)
+ if ( (lastDepth != -1) && (lastDepth > 0 || _bedGraphAll) ) {
cout << chrom << "\t" << lastStart << "\t" << pos << "\t" << lastDepth << endl;
}
//Set current position as the new interval start + depth
@@ -340,13 +336,12 @@ void BedGenomeCoverage::ReportChromCoverageBedGraph(const vector<DEPTH> &chromCo
lastStart = pos;
}
// Default: the depth has not changed, so we will not print anything.
- // Proceed until the depth changes.
-
+ // Proceed until the depth changes.
// Update depth
depth = depth - chromCov[pos].ends;
}
//Print information about the last position
- if (lastDepth != -1) {
+ if ( (lastDepth != -1) && (lastDepth > 0 || _bedGraphAll) ) {
cout << chrom << "\t" << lastStart << "\t" << chromSize << "\t" << lastDepth << endl;
}
}
diff --git a/src/genomeCoverageBed/genomeCoverageBed.h b/src/genomeCoverageBed/genomeCoverageBed.h
index 8b36de0d..a7ab6e5f 100755
--- a/src/genomeCoverageBed/genomeCoverageBed.h
+++ b/src/genomeCoverageBed/genomeCoverageBed.h
@@ -40,7 +40,7 @@ public:
// constructor
BedGenomeCoverage(string bedFile, string genomeFile, bool eachBase, bool startSites,
- bool bedGraph, int max, bool bamInput);
+ bool bedGraph, bool bedGraphAll, int max, bool bamInput);
// destructor
~BedGenomeCoverage(void);
@@ -54,6 +54,7 @@ private:
bool _eachBase;
bool _startSites;
bool _bedGraph;
+ bool _bedGraphAll;
int _max;
// The BED file from which to compute coverage.
diff --git a/src/genomeCoverageBed/genomeCoverageMain.cpp b/src/genomeCoverageBed/genomeCoverageMain.cpp
index 3dcf152b..3c725227 100755
--- a/src/genomeCoverageBed/genomeCoverageMain.cpp
+++ b/src/genomeCoverageBed/genomeCoverageMain.cpp
@@ -34,12 +34,13 @@ int main(int argc, char* argv[]) {
string genomeFile;
int max = 999999999;
- bool haveBed = false;
- bool bamInput = false;
- bool haveGenome = false;
- bool startSites = false;
- bool bedGraph = false;
- bool eachBase = false;
+ bool haveBed = false;
+ bool bamInput = false;
+ bool haveGenome = false;
+ bool startSites = false;
+ bool bedGraph = false;
+ bool bedGraphAll = false;
+ bool eachBase = false;
// check to see if we should print out some help
if(argc <= 1) showHelp = true;
@@ -87,7 +88,10 @@ int main(int argc, char* argv[]) {
}
else if(PARAMETER_CHECK("-bg", 3, parameterLength)) {
bedGraph = true;
- }
+ }
+ else if(PARAMETER_CHECK("-bga", 4, parameterLength)) {
+ bedGraphAll = true;
+ }
else if(PARAMETER_CHECK("-max", 4, parameterLength)) {
if ((i+1) < argc) {
max = atoi(argv[i + 1]);
@@ -106,14 +110,19 @@ int main(int argc, char* argv[]) {
showHelp = true;
}
if (bedGraph && eachBase) {
- cerr << endl << "*****" << endl << "*****ERROR: Use -d or -bedgraph, not both" << endl << "*****" << endl;
+ cerr << endl << "*****" << endl << "*****ERROR: Use -d or -bg, not both" << endl << "*****" << endl;
showHelp = true;
}
-
+ if (bedGraphAll && eachBase) {
+ cerr << endl << "*****" << endl << "*****ERROR: Use -d or -bga, not both" << endl << "*****" << endl;
+ showHelp = true;
+ }
+
if (!showHelp) {
BedGenomeCoverage *bc = new BedGenomeCoverage(bedFile, genomeFile, eachBase,
- startSites, bedGraph, max, bamInput);
+ startSites, bedGraph, bedGraphAll,
+ max, bamInput);
delete bc;
return 0;
@@ -127,9 +136,11 @@ void ShowHelp(void) {
cerr << endl << "Program: " << PROGRAM_NAME << " (v" << VERSION << ")" << endl;
- cerr << "Author: Aaron Quinlan (aaronquinlan@gmail.com)" << endl;
+ cerr << "Authors: Aaron Quinlan (aaronquinlan@gmail.com)" << endl;
+ cerr << " Gordon Assaf, CSHL" << endl;
+ cerr << " Royden Clark, UVa." << endl << endl;
- cerr << "Summary: Compute the coverage of a BED file among a genome." << endl << endl;
+ cerr << "Summary: Compute the coverage of a BED/BAM file among a genome." << endl << endl;
cerr << "Usage: " << PROGRAM_NAME << " [OPTIONS] -i <bed> -g <genome>" << endl << endl;
@@ -144,6 +155,10 @@ void ShowHelp(void) {
cerr << "\t-bg\t" << "Report depth in BedGraph format. For details, see:" << endl;
cerr << "\t\tgenome.ucsc.edu/goldenPath/help/bedgraph.html" << endl << endl;
+ cerr << "\t-bga\t" << "Report depth in BedGraph format (same as above)." << endl;
+ cerr << "\t\tHowever with this option, regions with zero " << endl;
+ cerr << "\t\tcoverage are also reported." << endl << endl;
+
cerr << "\t-max\t" << "Combine all positions with a depth >= max into" << endl;
cerr << "\t\ta single bin in the histogram. Irrelevant" << endl;
cerr << "\t\tfor -d and -bedGraph" << endl;
--
GitLab