Bug: HMMRATAC --save-digested outputs four identical signal tracks
mrendleman opened this issue · 0 comments
mrendleman commented
Describe the bug
When running MACS3 in HMMRATAC mode, using the --save-digested flag writes the short signal to all four signal tracks.
To Reproduce
Use the --save-digested flag and then run diff
on the resulting files.
Additional context
The problem is here:
https://github.com/macs3-project/MACS/blob/429f8f8b0114974b168900ba17d96467c27b550f/MACS3/Commands/hmmratac_cmd.py#L158C1-L170C21
Solution
The index of digested_atac_signals
was mistakenly set to 0 in all four cases. The following block:
if options.save_digested:
fhd = open(short_bdgfile,"w")
digested_atac_signals[ 0 ].write_bedGraph(fhd, "short","short")
fhd.close()
fhd = open(mono_bdgfile,"w")
digested_atac_signals[ 0 ].write_bedGraph(fhd, "mono","mono")
fhd.close()
fhd = open(di_bdgfile,"w")
digested_atac_signals[ 0 ].write_bedGraph(fhd, "di","di")
fhd.close()
fhd = open(tri_bdgfile,"w")
digested_atac_signals[ 0 ].write_bedGraph(fhd, "tri","tri")
fhd.close()
should actually be:
if options.save_digested:
fhd = open(short_bdgfile,"w")
digested_atac_signals[ 0 ].write_bedGraph(fhd, "short","short")
fhd.close()
fhd = open(mono_bdgfile,"w")
digested_atac_signals[ 1 ].write_bedGraph(fhd, "mono","mono")
fhd.close()
fhd = open(di_bdgfile,"w")
digested_atac_signals[ 2 ].write_bedGraph(fhd, "di","di")
fhd.close()
fhd = open(tri_bdgfile,"w")
digested_atac_signals[ 3 ].write_bedGraph(fhd, "tri","tri")
fhd.close()