macs3-project/MACS

Bug: HMMRATAC --save-digested outputs four identical signal tracks

mrendleman opened this issue · 0 comments

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()