From f66261fbdbef3b0683f8d8ef37a18f2d0baa554c Mon Sep 17 00:00:00 2001 From: Brent Pedersen Date: Sun, 4 Dec 2022 11:31:17 -0800 Subject: [PATCH] fix coverage formatting when scale == 1 closes #371 --- src/genomeCoverageBed/genomeCoverageBed.cpp | 27 +++++++++++++++++--- src/genomeCoverageBed/genomeCoverageMain.cpp | 6 ++--- test/genomecov/mk-deep.py | 7 +++++ test/genomecov/test-genomecov.sh | 7 +++++ 4 files changed, 41 insertions(+), 6 deletions(-) create mode 100644 test/genomecov/mk-deep.py diff --git a/src/genomeCoverageBed/genomeCoverageBed.cpp b/src/genomeCoverageBed/genomeCoverageBed.cpp index c6dafd732..77c0a7355 100644 --- a/src/genomeCoverageBed/genomeCoverageBed.cpp +++ b/src/genomeCoverageBed/genomeCoverageBed.cpp @@ -378,6 +378,12 @@ void BedGenomeCoverage::CoverageBam(string bamFile) { void BedGenomeCoverage::ReportChromCoverage(const vector &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 &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 &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 &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 &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 index 4a187975d..e94c3af52 100644 --- 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 -g OR -ibam " << 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 new file mode 100644 index 000000000..1e090fc7f --- /dev/null +++ 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 index 76df30ee7..5a269f211 100755 --- 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;