ac.tl.nucleosome_signal bug

[ad7115]

Describe the bug
I am running the muon tutorials for Single-cell RNA-seq and ATAC-seq integration on 10k PBMCs.
After running atac.tl.nucleosome_signal on the tutorial dataset, when plotting the histogram, I only get a square histogram, with a nucleosome_signal of 1.0 for every fragment. I tried the command without the barcode option (ac.tl.nucleosome_signal(atac, n=1e6)) and with it (ac.tl.nucleosome_signal(atac, n=1e6, barcodes="orig_barcode")) but nothing changed.

To Reproduce
I basically reproduced every step of the muon tutorial on this link : https://muon-tutorials.readthedocs.io/en/latest/single-cell-rna-atac/pbmc10k/1-Gene-Expression-Processing.html, until the ATAC-seq specific QC. Plot the histogramme with the command: mu.pl.histogram(atac, "nucleosome_signal", kde=False).

Expected behaviour
The expected histogram shown in the tutorial has a range of values between 0 and 2.5 approximately.

System

• Python 3.8

[russellgould]

I also ran into this issue. I did some digging and it turns out that the problem is stemming from this line combined with the try...except block catching ALL exceptions. In this case it's catching the informative AttributeError that tells us that the object returned by pysam.TabixFile().fetch() doesn't have a .next() function.

So, the code is just skipping every single fragment in the file since every call to f.next() raises an exception, which is then caught, and the iterator passes to the next iteration.

The simplest way to fix this function is to copy the function to your code and just change this line:

f = fr.next()

to this:

f = next(fr)

A safer solution would be to change the try...except block to only catch the KeyError exceptions for fragments that don't exist in the mdata object. This would allow other errors to be raised when things are actually broken.

for _ in tqdm(range(n), desc="Reading Fragments"):
    try:
        f = next(fr) # changed from fr.next()
        length = f.end - f.start
        row_ind = d[f.name]
        if length < nucleosome_free_upper_bound:
            mat[row_ind, 0] += 1
        elif length < mononuleosomal_upper_bound:
            mat[row_ind, 1] += 1
    except KeyError: # changed to only catch KeyError (safer for debugging)
        pass

[gtca]

Thanks a lot for noticing and debugging this, @ad7115 and @russellgould!

I've updated the function in the suggested changes, and it should be out in v0.1.4. I am curious if this might have been something Windows-specific as I wasn't able to reproduce it on Unix...

[russellgould]

Great! Thanks for making that change!

I'm on a Mac so should also give an error on other flavors of Unix. Possibly we have a different version of pysam installed? I found the above error running muon version 0.1.3 and pysam version 0.21.0.