maxatac prepare crashes on chrX because ATAC_bowtie2_pipeline.sh ignores the --chromosomes argument

[PanosFirmpas]

When running maxatac prepare on a bulk ATAC-seq BAM file with --chromosomes including chrX, the run crashes during normalization with:

RuntimeError: The end coordinate must be a number!

The root cause is in ATAC_bowtie2_pipeline.sh (in MiraldiLab/maxATAC_data, at scripts/ATAC/ATAC_bowtie2_pipeline.sh), not in the Python normalization code where the traceback points.

Root cause

maxatac/analyses/prepare.py passes the user's chromosome list as an extra positional argument to ATAC_bowtie2_pipeline.sh:

subprocess.run(["bash", PREPARE_BULK_SCRIPT, args.input, args.name, output_dir,
                 str(args.threads), args.blacklist, args.chrom_sizes, str(args.slop),
                 str(scale_factor), "deduplicate", " ".join(args.chromosomes)], check=True)

But ATAC_bowtie2_pipeline.sh only declares and reads parameters $1 through $9 — it never reads the 11th argument (the chromosome list). Instead it uses a hardcoded, autosome-only variable:

keepChr='chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr20 chr21 chr22'

which is used to filter reads via samtools view:

samtools view -@ ${cores} -f 3 -b ${deduped} ${keepChr} | ...

So even when a user passes --chromosomes ... chrX, chrX reads are silently dropped before the BAM is converted to a bedgraph/bigwig. The resulting .bw file ends up with no chrX contig in its header at all.

Later, get_genomic_stats() in maxatac/utilities/normalization_tools.py iterates over the user-requested chromosome list (built independently from the .chrom.sizes file) and calls:

chr_vals = np.nan_to_num(input_bigwig.values(chromosome, 0, input_bigwig.chroms(chromosome), numpy=True))

Since chrX isn't present in the bigwig's header, pyBigWig's chroms("chrX") returns None (this is documented pyBigWig behavior for absent chromosomes), which is then passed as the end coordinate to .values(), producing the crash.