Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fix coverage formatting when scale == 1 #1031

Merged
merged 1 commit into from
Dec 6, 2022
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
27 changes: 24 additions & 3 deletions src/genomeCoverageBed/genomeCoverageBed.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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++) {
Expand All @@ -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);
Expand Down Expand Up @@ -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;
Expand All @@ -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
Expand All @@ -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);
}
}
6 changes: 3 additions & 3 deletions src/genomeCoverageBed/genomeCoverageMain.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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;
Expand All @@ -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;

Expand Down
7 changes: 7 additions & 0 deletions test/genomecov/mk-deep.py
Original file line number Diff line number Diff line change
@@ -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*')
7 changes: 7 additions & 0 deletions test/genomecov/test-genomecov.sh
Original file line number Diff line number Diff line change
Expand Up @@ -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;