Runtime question - atlas-scale analysis

Hello! I love this package, so thank you so much for a python-native solution to scATAC analysis. It's been indispensable.

I had a question about runtime for the atlas-level dataset. I was reading through the preprint and noted this:

This atlas comprises 92 scATAC-seq139
samples, roughly 650,000 cells, over 23 billion raw reads, and has a total size of 1.6 terabytes. As shown in Fig. 1e,140
SnapATAC2 completed the entire process in 5.22 hours, whereas ArchR took 14.27 hours. At this scale, SnapATAC2141
is almost three times faster than ArchR, decreasing computational costs by approximately 63.4%

I wanted to actually use this data for some analysis, and so I grabbed it from geo, which was linked in the original paper for that dataset. I've also got the metadata.

Using the tutorial for analyzing this dataset, I developed this code:

import os
import multiprocessing as mp
from pathlib import Path
from glob import glob

import snapatac2 as snap
import numpy as np

from tqdm import tqdm

snap.__version__

data_dir = os.path.expandvars("$DATA/geniml/datasets/scatlas/fragments")
output_dir = os.path.expandvars("$DATA/geniml/datasets/scatlas/h5ad")

if not os.path.exists(output_dir):
    os.makedirs(output_dir)

fragment_files = glob(os.path.join(data_dir, "*.bed.gz"))

outputs = []
for file in fragment_files:
    new_file = Path(file).name.replace(".bed.gz", ".h5ad")
    output_file = os.path.join(output_dir, new_file)
    outputs.append(output_file)

adatas = snap.pp.import_data(
    fragment_files, 
    file=outputs, 
    chrom_sizes=snap.genome.hg38,
    min_num_fragments=1000, 
)

snap.pp.add_tile_matrix(adatas, bin_size=5000, n_jobs=mp.cpu_count() - 2)
snap.pp.select_features(adatas, n_jobs=mp.cpu_count() - 2)
snap.pp.scrublet(adatas, n_jobs=mp.cpu_count() - 2)
snap.pp.filter_doublets(adatas, n_jobs=mp.cpu_count() - 2)

adataset = snap.AnnDataSet(
    adatas=[(f.filename.split('/')[-1].split('.h5ad')[0], f) for f in adatas],
    filename="data.h5ads"
)

adataset.obs_names = np.array(adataset.obs['sample']) + '+' + np.array(adataset.obs_names)

snap.pp.select_features(adataset, n_features=250000)
snap.tl.spectral(adataset)

# Store tissue types in .obs
adataset.obs['tissue'] = [x.split(':')[0] for x in adataset.obs['sample']]
snap.pp.mnc_correct(adataset, batch="sample", groupby='tissue', key_added='X_spectral')

Everything runs nicely, but the snap.pp.scrublet snippet predicts it won't finish for another 20 hours, which is quite long. Was wondering if there was something I was doing wrong - how important is the n_jobs param? I'm working with 40 cores + 64G RAM.

Any help is appreciated! Thank you.

Start with the atlas tutorial and the snap.pp.scrublet call in the issue, then inspect how its n_jobs parameter uses the available 40 cores and 64G RAM. Compare the reported runtime with the surrounding preprocessing steps and document whether the invocation is expected; done means explaining the runtime or identifying a reproducible configuration problem.