The maxatac prepare function was initially created as a convenience function for filtering, inferring Tn5 sites, and converting cut -site level coverage to min-max normalized bigwig tracks in one step. This function calls on the bash scripts that were used by our snakemake/cwl/bash workflows for ATAC-seq data processing. This will most likely be how most users prepare data as opposed to going through each step individually, so we should think about improving the user experience.

[x] The biggest bug we should work on is adding error reporting and some code to catch if the proper outputs were not produced when using python to run our shell scripts. @anthonybejjani was using maxatac prepare to prepare scATAC-seq fragment files. He was able to run the script to completion and did not get an error message that: 1) A problem with bedGraphToBigWig finding the correct libssl library #102
2) All of the expected output files were not produced.

The problem is that it appears to the user that the run completes correctly, despite having encountered some error during running the shell script. We should add more logging information during processing. We should also look into whether we should use python to execute the shell commands as opposed to just running a shell script from python with the commands internally. We could also add code to the shell script to catch problems with execution.

[x] We should also add a test for maxatac functions that will make sure pybigwig is installed and can find numpy correctly, before running through the entire workflow. This is related to #96. We should point to the fix if the issue is detected.
[x] Double check that all of the unnecessary bedgraphs and intermediate files are removed to save space. We might want to add flags for whether to save specific intermediate files.
[x] Update and add better logging messages for different processes running. At least add messages for major events like saving files or removing files. We could also have a final printout that has the names and locations of all files and their sizes.

Hey, I gave maxatac prepare this input:

`#!/bin/bash

BSUB -W 12:00

BSUB -M 200000

BSUB -n 24

BSUB -J kessi_maxATAC_prepare

BSUB -o /data/miraldiLab/team/kessi/out_files/maxATAC_prepare_test_071522.out

BSUB -e /data/miraldiLab/team/kessi/err_files/maxATAC_prepare_test_071522.err

module load anaconda

module load anaconda (I don't remember the exact number but I think it is 3 you can check with module avail)

module load anaconda3/1.0.0

conda activate your env

conda activate YOUR_CONDA_ENV_NAME # This env should have maxatac installed

source activate maxatac_test

No need for the sh infront of maxatac

maxatac prepare \ -i /data/miraldiLab/team/kessi/snakeATAC/output_dir/Th0_48h_r1/alignments/Th0_48h_r1.bam \ -o /data/miraldiLab/team/kessi/maxATAC \ -prefix Th0_48h_r1 \ --blacklist_bed /users/dia6sx/opt/maxatac/data/mm10/mm10_maxatac_blacklist.bed \ --blacklist_bw /users/dia6sx/opt/maxatac/data/mm10/mm10_maxatac_blacklist.bw \ --chrom_sizes /users/dia6sx/opt/maxatac/data/mm10/mm10.chrom.sizes \ -chroms chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19`

`[2022-07-15 16:28:03,370] Input file: /data/miraldiLab/team/kessi/snakeATAC/output_dir/Th0_48h_r1/alignments/Th0_48h_r1.bam Input chromosome sizes file: /users/dia6sx/opt/maxatac/data/mm10/mm10.chrom.sizes Tn5 cut sites will be slopped 20 bps on each side Input blacklist file: /users/dia6sx/opt/maxatac/data/hg38/hg38_maxatac_blacklist.bw Output filename: Th0_48h_r1 Output directory: /data/miraldiLab/team/kessi/maxATAC Using a millions factor of: 20000000 Using 12 threads to run job. [2022-07-15 16:28:03,471] Generate the normalized signal tracks. [2022-07-15 16:28:03,471] Working on a bulk ATAC-seq BAM file Getting the number of reads in the BAM file [2022-07-15 16:28:10,236] Processing BAM to bigwig. Running eduplication fixmate: invalid option -- '@' Usage: samtools fixmate Options: -r Remove unmapped reads and secondary alignments -p Disable FR proper pair check -c Add template cigar ct tag --input-fmt-option OPT[=VAL] Specify a single input file format option in the form of OPTION or OPTION=VALUE -O, --output-fmt FORMAT[,OPT[=VAL]]... Specify output format (SAM, BAM, CRAM) --output-fmt-option OPT[=VAL] Specify a single output file format option in the form of OPTION or OPTION=VALUE --reference FILE Reference sequence FASTA FILE [null]

As elsewhere in samtools, use '-' as the filename for stdin/stdout. The input file must be grouped by read name (e.g. sorted by name). Coordinated sorted input is not accepted. index: invalid option -- '@' Usage: samtools index [-bc] [-m INT] [out.index] Options: -b Generate BAI-format index for BAM files [default] -c Generate CSI-format index for BAM files -m INT Set minimum interval size for CSI indices to 2^INT [14] [main] unrecognized command 'markdup' index: invalid option -- '@' Usage: samtools index [-bc] [-m INT] [out.index] Options: -b Generate BAI-format index for BAM files [default] -c Generate CSI-format index for BAM files -m INT Set minimum interval size for CSI indices to 2^INT [14] rm: cannot remove ‘Th0_48h_r1_fixmate.bam.bai’: No such file or directory [main_samview] random alignment retrieval only works for indexed BAM or CRAM files. index: invalid option -- '@' Usage: samtools index [-bc] [-m INT] [out.index] Options: -b Generate BAI-format index for BAM files [default] -c Generate CSI-format index for BAM files -m INT Set minimum interval size for CSI indices to 2^INT [14] rm: cannot remove ‘Th0_48h_r1_deduped.bam.bai’: No such file or directory bedGraphToBigWig: error while loading shared libraries: libssl.so.1.0.0: cannot open shared object file: No such file or directory [2022-07-15 16:29:36,562] Min-max normalize signal tracks [2022-07-15 16:29:36,565] Normalization Input bigwig file: /data/miraldiLab/team/kessi/maxATAC/Th0_48h_r1_IS_slop20_RP20M.bw Output filename: /data/miraldiLab/team/kessi/maxATAC/Th0_48h_r1_IS_slop20_RP20M_minmax01.bw Output directory: /data/miraldiLab/team/kessi/maxATAC Using min-max normalization [2022-07-15 16:29:36,583] Calculating stats per chromosome [urlOpen] Couldn't open /data/miraldiLab/team/kessi/maxATAC/Th0_48h_r1_IS_slop20_RP20M.bw for reading [urlOpen] Couldn't open /data/miraldiLab/team/kessi/maxATAC/Th0_48h_r1_IS_slop20_RP20M.bw for reading [pyBwOpen] bw is NULL! Traceback (most recent call last): File "/users/dia6sx/.conda/envs/maxatac_test/bin/maxatac", line 24, in sys.exit(main(sys.argv[1:])) File "/users/dia6sx/.conda/envs/maxatac_test/bin/maxatac", line 20, in main args.func(args) File "/users/dia6sx/.conda/envs/maxatac_test/lib/python3.9/site-packages/maxatac/analyses/prepare.py", line 142, in run_prepare run_normalization(args) File "/users/dia6sx/.conda/envs/maxatac_test/lib/python3.9/site-packages/maxatac/analyses/normalize.py", line 64, in run_normalization min_value, max_value, median, mad, mean_value, std_value = get_genomic_stats(bigwig_path=args.signal, File "/users/dia6sx/.conda/envs/maxatac_test/lib/python3.9/site-packages/maxatac/utilities/normalization_tools.py", line 24, in get_genomic_stats with pyBigWig.open(bigwig_path) as input_bigwig: RuntimeError: Received an error during file opening!`

There seems to be a problem when preparing a bam file from snakeATAC. Note I am able to bypass the prepare step by using the bw generated.

MiraldiLab / maxATAC

Improve `maxatac prepare` #107