fix coverage formatting when scale == 1

closes #371
arq5x · Dec 4, 2022 · f66261f · f66261f
1 parent 4c19da5
commit f66261f
Show file tree

Hide file tree

Showing 4 changed files with 41 additions and 6 deletions.
diff --git a/src/genomeCoverageBed/genomeCoverageBed.cpp b/src/genomeCoverageBed/genomeCoverageBed.cpp
@@ -378,6 +378,12 @@ void BedGenomeCoverage::CoverageBam(string bamFile) {
 void BedGenomeCoverage::ReportChromCoverage(const vector<DEPTH> &chromCov, const CHRPOS &chromSize, const string &chrom, chromHistMap &chromDepthHist) {
 
     if (_eachBase) {
+        streamsize prec;
+        std::ios_base::fmtflags fflags;
+        if(_scale == 1.0) {
+            prec = cout.precision(0);
+            fflags = cout.setf(std::ios_base::fixed);
+        }
         int depth = 0; // initialize the depth
         CHRPOS offset = (_eachBaseZeroBased)?0:1;
         for (CHRPOS pos = 0; pos < chromSize; pos++) {
@@ -386,11 +392,15 @@ void BedGenomeCoverage::ReportChromCoverage(const vector<DEPTH> &chromCov, const
             // report the depth for this position.
             if (depth>0 || !_eachBaseZeroBased) {
                 cout << chrom << "\t" << pos+offset << "\t"
-                     << (_scale == 1.0 ? depth : depth * _scale)
+                     << depth * _scale
                      << endl;
             }
             depth = depth - chromCov[pos].ends;
         }
+        if(_scale == 1.0) {
+            cout.precision(prec);
+            cout.setf(fflags);
+        }
     }
     else if (_bedGraph == true || _bedGraphAll == true) {
         ReportChromCoverageBedGraph(chromCov, chromSize, chrom);
@@ -472,11 +482,18 @@ void BedGenomeCoverage::ReportGenomeCoverage(chromHistMap &chromDepthHist) {
 }
 
 
+
 void BedGenomeCoverage::ReportChromCoverageBedGraph(const vector<DEPTH> &chromCov, const CHRPOS &chromSize, const string &chrom) {
 
     int depth = 0; // initialize the depth
     CHRPOS lastStart = -1;
     int lastDepth = -1;
+    streamsize prec;
+    std::ios_base::fmtflags fflags;
+    if(_scale == 1.0) {
+        prec = cout.precision(0);
+        fflags = cout.setf(std::ios_base::fixed);
+    }
 
     for (CHRPOS pos = 0; pos < chromSize; pos++) {
         depth += chromCov[pos].starts;
@@ -497,7 +514,7 @@ void BedGenomeCoverage::ReportChromCoverageBedGraph(const vector<DEPTH> &chromCo
                      << "\t"
                      << pos
                      << "\t"
-                     << (_scale == 1.0 ? lastDepth : lastDepth * _scale)
+                     << lastDepth * _scale
                      << endl;
             }
             //Set current position as the new interval start + depth
@@ -517,7 +534,11 @@ void BedGenomeCoverage::ReportChromCoverageBedGraph(const vector<DEPTH> &chromCo
              << "\t"
              << chromSize
              << "\t"
-             << (_scale == 1.0 ? lastDepth : lastDepth * _scale)
+             << lastDepth * _scale
              << endl;
     }
+    if(_scale == 1.0) {
+        cout.precision(prec);
+        cout.setf(fflags);
+    }
 }
diff --git a/src/genomeCoverageBed/genomeCoverageMain.cpp b/src/genomeCoverageBed/genomeCoverageMain.cpp
@@ -220,7 +220,7 @@ int genomecoverage_main(int argc, char* argv[]) {
       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,
@@ -241,7 +241,7 @@ int genomecoverage_main(int argc, char* argv[]) {
 void genomecoverage_help(void) {
 
     cerr << "\nTool:    bedtools genomecov (aka genomeCoverageBed)" << endl;
-    cerr << "Version: " << VERSION << "\n";    
+    cerr << "Version: " << VERSION << "\n";
     cerr << "Summary: Compute the coverage of a feature file among a genome." << endl << endl;
 
     cerr << "Usage: " << PROGRAM_NAME << " [OPTIONS] -i <bed/gff/vcf> -g <genome> OR -ibam <bam/cram>" << endl << endl;
@@ -253,7 +253,7 @@ void genomecoverage_help(void) {
 
     cerr << "\t-g\t\t" << "Provide a genome file to define chromosome lengths." << endl;
     cerr << "\t\t\tNote:Required when not using -ibam option." << endl << endl;
-    
+
     cerr << "\t-d\t\t" << "Report the depth at each genome position (with one-based coordinates)." << endl;
     cerr << "\t\t\tDefault behavior is to report a histogram." << endl << endl;
 

diff --git a/test/genomecov/mk-deep.py b/test/genomecov/mk-deep.py
@@ -0,0 +1,7 @@
+print("""\
+@HD\tVN:1.0\tSO:coordinate
+@SQ\tSN:c1\tAS:genome.txt\tLN:100""")
+
+#y1      0       1       16      255     5M      *       0       0       *       *
+for i in range(1000000):
+    print(f'r{i}\t0\tc1\t1\t100\t100M\t*\t0\t0\t*\t*')
diff --git a/test/genomecov/test-genomecov.sh b/test/genomecov/test-genomecov.sh
@@ -288,4 +288,11 @@ CRAM_REFERENCE=test_ref.fa $BT genomecov -ibam empty.cram > obs
 check obs exp
 rm obs exp
 
+python mk-deep.py > deep.sam
+echo -e "    genomecov.t18...\c"
+echo "c1	1	1000000" > exp
+$BT genomecov -d -ibam deep.sam | head -1 > obs
+check obs exp
+rm obs exp deep.sam
+
 [[ $FAILURES -eq 0 ]] || exit 1;