macs2 error, it worked fine before

[danielcgingerich]

macs2 in signac has worked fine before. Now I am receiving an error "ZeroDivisionError: float division by zero". I am using the parameters idents = and group.by = to call peaks for one cell type at a time. From what I read on biostars, my guess is that the new fragment file that was created is missing an index. But why is this happening? The code I am running was not changed, and it worked fine just yesterday.

library(Signac)
library(Seurat)


#project skeleton\~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
number.list<- c("99", "111", "127", "191", "196",
                "273", "347", "357", "372", "430",
                "542", "601","676", "688", "730", 
                "781", "963", "984", "1099", "1545",
                "1557", "1600", "1670", "1690")
main.dir <- '/gpfs/fs1/data/chiba-faleklab/Daniel/snATAC.snRNA.R01.v2/R.objects/'
rds.list <- paste0('rds', c(97:120), '_', number.list, 'srt_preprocessed_filtered.rds')
meta.list <- paste0('rds', c(249:272), '_', number.list, 'predicted.id.celltype_snRNA_bugfix.rds')
endos.vlmcs <- paste0('rds', c(173:196), '_', number.list, 'predicted.id_bugfix.rds')
celltype <- c('Exc', 'Inh', 'Astro', 'Micro', 'Oligo', 'OPC')
celltype2 <- c('exc2', 'inh2', 'astro2', 'micro2', 'oligo2', 'opc2')

slurm.task.id <- Sys.getenv('SLURM_ARRAY_TASK_ID')
i <- as.numeric(slurm.task.id)

j = 1


#load objects\~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
outs <- paste0(main.dir, 'recalled.peaks/', celltype2[j], '/', number.list[i], '/')
atac <- readRDS(paste0(main.dir, rds.list[i]))
metadata <- readRDS(paste0(main.dir, meta.list[i]))
endos.vlmcs <- readRDS(paste0(main.dir, endos.vlmcs[i]))
atac <- AddMetaData(atac, metadata = metadata)

#call peaks\~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
message('\ncalling peaks \n')
message(paste0('dataset: ', number.list[i]))
peaks <- CallPeaks(object = atac, 
                   macs2.path = 'macs2', 
                   idents = celltype[j],
                   group.by = 'predicted.id', 
                   combine.peaks = FALSE,
                   cleanup = FALSE,
                   outdir = outs)
Processing file /gpfs/fs1/data/chiba-faleklab/snATAC_fastq_count_order6394_07232020/HMWT7DRXX/99atac_count/outs/fragments.tsv.gz
Keeping 574 cell barcodes
Splitting into 1 files
Done Processing 259 million lines

INFO  @ Tue, 20 Apr 2021 15:14:51: 
# Command line: callpeak -t /tmp/RtmptVilAk/Exc.bed -g 2.7e+09 -f BED --nomodel --extsize 200 --shift -100 -n Exc --outdir /gpfs/fs1/data/chiba-faleklab/Daniel/snATAC.snRNA.R01.v2/R.objects/recalled.peaks/exc2/99/
# ARGUMENTS LIST:
# name = Exc
# format = BED
# ChIP-seq file = ['/tmp/RtmptVilAk/Exc.bed']
# control file = None
# effective genome size = 2.70e+09
# band width = 300
# model fold = [5, 50]
# qvalue cutoff = 5.00e-02
# The maximum gap between significant sites is assigned as the read length/tag size.
# The minimum length of peaks is assigned as the predicted fragment length "d".
# Larger dataset will be scaled towards smaller dataset.
# Range for calculating regional lambda is: 10000 bps
# Broad region calling is off
# Paired-End mode is off
 
INFO  @ Tue, 20 Apr 2021 15:14:51: #1 read tag files... 
INFO  @ Tue, 20 Apr 2021 15:14:51: #1 read treatment tags... 
Exception ZeroDivisionError: 'integer division or modulo by zero' in 'MACS2.IO.Parser.GenericParser.tsize' ignored
INFO  @ Tue, 20 Apr 2021 15:14:51: #1 tag size is determined as 0 bps 
INFO  @ Tue, 20 Apr 2021 15:14:51: #1 tag size = 0.0 
INFO  @ Tue, 20 Apr 2021 15:14:51: #1  total tags in treatment: 0 
INFO  @ Tue, 20 Apr 2021 15:14:51: #1 user defined the maximum tags... 
INFO  @ Tue, 20 Apr 2021 15:14:51: #1 filter out redundant tags at the same location and the same strand by allowing at most 1 tag(s) 
INFO  @ Tue, 20 Apr 2021 15:14:51: #1  tags after filtering in treatment: 0 
Traceback (most recent call last):
  File "/nfs/software/helmod/apps/Core/macs2/2.1.2.20181108-gcb01/bin/macs2", line 622, in <module>
    main()
  File "/nfs/software/helmod/apps/Core/macs2/2.1.2.20181108-gcb01/bin/macs2", line 57, in main
    run( args )
  File "/nfs/software/helmod/apps/Core/macs2/2.1.2.20181108-gcb01/lib/python2.7/site-packages/MACS2/callpeak_cmd.py", line 112, in run
    info("#1  Redundant rate of treatment: %.2f", float(t0 - t1) / t0)
ZeroDivisionError: float division by zero
Error in file(file, "rt") : cannot open the connection
Calls: CallPeaks ... CallPeaks -> CallPeaks.default -> read.table -> file
In addition: Warning message:
In file(file, "rt") :
  cannot open file '/gpfs/fs1/data/chiba-faleklab/Daniel/snATAC.snRNA.R01.v2/R.objects/recalled.peaks/exc2/99//Exc_peaks.narrowPeak': No such file or directory
Execution halted

[timoast]

@danielcgingerich the index is not actually needed for MACS2 peak calling. What changed between when the code worked and did not work? Did you update Signac, or update any other code involved here?

[danielcgingerich]

@timoast no I did not change anything or update signac. I was submitting this script as a slurm job array, hence the code:

slurm.task.id <- Sys.getenv('SLURM_ARRAY_TASK_ID')
i <- as.numeric(slurm.task.id)

There are 24 instances of this script running in parallel for the array. I believe the issue may be with the creation of so many tmp directories at once. My understanding is that R can only make one at a time. When the array is running, this would mean 24 tmp directories are created at once. Not sure if running R through a slurm array on my university hpc would allow for the creation of multiple tmp directories. Would you know?

I fixed the issue by running macs2 through R as written in your source code, but instead of using tmp directories, I created permanent directories ahead of time to store the .bed files.

[timoast]

Ok it sounds like this might be specific to the particular way you're running it on a cluster, but we could consider adding a parameter for the path to store the split fragment files that should also solve this issue.

[timoast]

Added a fragment.tempdir parameter to CallPeaks()