First cut at -scale for genomeCoverageBed

src/genomeCoverageBed/genomeCoverageBed.cpp

@@ -13,9 +13,11 @@ Licenced under the GNU General Public License 2.0 license.
 #include "genomeCoverageBed.h"
 
 
-BedGenomeCoverage::BedGenomeCoverage(string bedFile, string genomeFile, bool eachBase,
-                                     bool startSites, bool bedGraph, bool bedGraphAll,
-                                     int max, bool bamInput, bool obeySplits,
+BedGenomeCoverage::BedGenomeCoverage(string bedFile, string genomeFile, 
+                                     bool eachBase, bool startSites, 
+                                     bool bedGraph, bool bedGraphAll,
+                                     int max, float scale,
+                                     bool bamInput, bool obeySplits,
                                      bool filterByStrand, string requestedStrand,
                                      bool only_5p_end, bool only_3p_end,
                                      bool eachBaseZeroBased,
@@ -29,6 +31,7 @@ BedGenomeCoverage::BedGenomeCoverage(string bedFile, string genomeFile, bool eac
     _bedGraph = bedGraph;
     _bedGraphAll = bedGraphAll;
     _max = max;
+    _scale = scale;
     _bamInput = bamInput;
     _obeySplits = obeySplits;
     _filterByStrand = filterByStrand;
@@ -269,7 +272,7 @@ void BedGenomeCoverage::ReportChromCoverage(const vector<DEPTH> &chromCov, const
             depth += chromCov[pos].starts;
             // report the depth for this position.
             if (depth>0 || !_eachBaseZeroBased)
-                cout << chrom << "\t" << pos+offset << "\t" << depth << endl;
+                cout << chrom << "\t" << pos+offset << "\t" << depth * _scale << endl;
             depth = depth - chromCov[pos].ends;
         }
     }
@@ -368,7 +371,7 @@ void BedGenomeCoverage::ReportChromCoverageBedGraph(const vector<DEPTH> &chromCo
             // (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;
+                cout << chrom << "\t" << lastStart << "\t" << pos << "\t" << lastDepth * _scale << endl;
             }
             //Set current position as the new interval start + depth
             lastDepth = depth;
@@ -381,6 +384,6 @@ void BedGenomeCoverage::ReportChromCoverageBedGraph(const vector<DEPTH> &chromCo
     }
     //Print information about the last position
     if ( (lastDepth != -1) && (lastDepth > 0 || _bedGraphAll) ) {
-        cout << chrom << "\t" << lastStart << "\t" << chromSize << "\t" << lastDepth << endl;
+        cout << chrom << "\t" << lastStart << "\t" << chromSize << "\t" << lastDepth * _scale << endl;
     }
 }

src/genomeCoverageBed/genomeCoverageBed.h

@@ -41,8 +41,11 @@ class BedGenomeCoverage {
 public:
 
     // constructor
-    BedGenomeCoverage(string bedFile, string genomeFile, bool eachBase, bool startSites,
-                      bool bedGraph, bool bedGraphAll, int max, bool bamInput, bool obeySplits,
+    BedGenomeCoverage(string bedFile, string genomeFile, 
+                      bool eachBase, bool startSites,
+                      bool bedGraph, bool bedGraphAll, 
+                      int max, float scale, 
+                      bool bamInput, bool obeySplits,
                       bool filterByStrand, string requestedStrand,
                       bool only_5p_end, bool only_3p_end,
                       bool eachBaseZeroBased,
@@ -63,6 +66,7 @@ private:
     bool _bedGraph;
     bool _bedGraphAll;
     int _max;
+    float _scale;
     bool _obeySplits;
     bool _filterByStrand;
     bool _only_5p_end;

src/genomeCoverageBed/genomeCoverageMain.cpp

@@ -33,6 +33,7 @@ int main(int argc, char* argv[]) {
     string bedFile;
     string genomeFile;
     int max = INT_MAX;
+    float scale = 1.0;
 
     bool haveBed = false;
     bool bamInput = false;
@@ -43,6 +44,7 @@ int main(int argc, char* argv[]) {
     bool eachBase = false;
     bool eachBaseZeroBased = false;
     bool obeySplits = false;
+    bool haveScale = false;
     bool filterByStrand = false;
     bool only_5p_end = false;
     bool only_3p_end = false;
@@ -110,6 +112,13 @@ int main(int argc, char* argv[]) {
                 i++;
             }
         }
+        else if(PARAMETER_CHECK("-scale", 6, parameterLength)) {
+            if ((i+1) < argc) {
+                haveScale = true;
+                scale = atof(argv[i + 1]);
+                i++;
+            }
+        }
         else if(PARAMETER_CHECK("-split", 6, parameterLength)) {
             obeySplits = true;
         }
@@ -182,11 +191,16 @@ int main(int argc, char* argv[]) {
       showHelp = true;
     }
 
+    if (haveScale && !(bedGraph||bedGraphAll||eachBase)) {
+      cerr << endl << "*****" << endl << "*****ERROR: Using -scale requires bedGraph output (use -bg or -bga) or per base depth (-d)." << endl << "*****" << endl;
+      showHelp = true;
+    }
+    
     if (!showHelp) {
 
         BedGenomeCoverage *bc = new BedGenomeCoverage(bedFile, genomeFile, eachBase,
                                                       startSites, bedGraph, bedGraphAll,
-                                                      max, bamInput, obeySplits,
+                                                      max, scale, bamInput, obeySplits,
                                                       filterByStrand, requestedStrand,
                                                       only_5p_end, only_3p_end,
                                                       eachBaseZeroBased,
@@ -252,6 +266,12 @@ void ShowHelp(void) {
     cerr << "\t\t\tfor -d and -bedGraph" << endl;
     cerr << "\t\t\t- (INTEGER)" << endl << endl;
 
+    cerr << "\t-scale\t\t" << "Scale the coverage by a constant factor." << endl;
+    cerr << "\t\t\tEach coverage value is multiplied by this factor before being reported." << endl;
+    cerr << "\t\t\tUseful for normalizing coverage by, e.g., reads per million (RPM)." << endl;
+    cerr << "\t\t\t- Default is 1.0; i.e., unscaled." << endl;
+    cerr << "\t\t\t- (FLOAT)" << endl << endl;
+
     cerr << "\t-trackline\t" << "Adds a UCSC/Genome-Browser track line definition in the first line of the output." << endl;
     cerr <<"\t\t\t- See here for more details about track line definition:" << endl;
     cerr <<"\t\t\t      http://genome.ucsc.edu/goldenPath/help/bedgraph.html" << endl;