issues decoding >5 modifications from MM/ML tags

[lpryszcz]

Hi, first of all thanks for developing pysam, it's fantastic package! I'm using the latest version (0.23.3).

Lately, I encountered several problems while working with modifications written by remora infer.

Everything works fine if there is only 1 modification in the BAM, for example

MM:Z:T+76?,0,0,0;       ML:B:C,4,5,4

But when I store several modifications (codes in range 76-94), pysam fails with ValueError: chr() arg not in range(0x110000).
Interestingly, both samtools and IGV are handling the modifications from the same BAM file correctly.

I tried to narrow it down and noticed when storing the same modifications as single letter codes m, h, c, f, g, e, b, a, o, n,
those are written in BAM as

MM:Z:T+m?,0;T+h?,0;T+c?,0;T+f?,0;T+g?,0;T+e?,0;T+b?,0;T+a?,0;T+o?,0;T+n?,0;

while pysam decode them incorrectly as

('T', 0, 'm'), ('T', 0, 'h'), ('T', 0, 'c'), ('T', 0, 'f'), ('T', 0, 'g'), ('\x00', -1098023120, '\x01'), ('嶨', -1137636000, 483813872), ('碴', 0, 1120917344), ('\x00', -1154113870, 0), ('碴', 0, 1099221952)

It seems, the issue is after 5th modification. I prepared two BAM files containing 6 modifications

MM:Z:T+m?,0;T+h?,0;T+f?,0;T+c?,0;T+g?,0;T+e?,0; ML:B:C,4,1,0,1,2,2
MM:Z:T+m?,0;T+h?,0;T+f?,0;T+c?,0;T+g?,0;T+e?,0; ML:B:C,3,0,0,0,0,0

and 5 modifications (g was skipped)

MM:Z:T+m?,0;T+h?,0;T+f?,0;T+c?,0;T+e?,0;        ML:B:C,4,1,0,1,2
MM:Z:T+m?,0;T+h?,0;T+f?,0;T+c?,0;T+e?,0;        ML:B:C,3,0,0,0,0

and pysam failed with 6 modifications

# print(a.get_tag('MM'), a.get_tag('ML'), a.modified_bases.keys())
T+m?,0;T+h?,0;T+f?,0;T+c?,0;T+g?,0;T+e?,0; array('B', [4, 1, 0, 1, 2, 2]) dict_keys([('T', 0, 'm'), ('T', 0, 'h'), ('T', 0, 'f'), ('T', 0, 'c'), ('T', 0, 'g'), ('\x00', -144053840, '\x1e')])
T+m?,0;T+h?,0;T+f?,0;T+c?,0;T+g?,0;T+e?,0; array('B', [3, 0, 0, 0, 0, 0]) dict_keys([('T', 0, 'm'), ('T', 0, 'h'), ('T', 0, 'f'), ('T', 0, 'c'), ('T', 0, 'g'), ('瑄', -144053840, 187320624)])

but processed correctly 5 modifications

# print(a.get_tag('MM'), a.get_tag('ML'), a.modified_bases.keys())
T+m?,0;T+h?,0;T+f?,0;T+c?,0;T+e?,0; array('B', [4, 1, 0, 1, 2]) dict_keys([('T', 0, 'm'), ('T', 0, 'h'), ('T', 0, 'f'), ('T', 0, 'c'), ('T', 0, 'e')])
T+m?,0;T+h?,0;T+f?,0;T+c?,0;T+e?,0; array('B', [3, 0, 0, 0, 0]) dict_keys([('T', 0, 'm'), ('T', 0, 'h'), ('T', 0, 'f'), ('T', 0, 'c'), ('T', 0, 'e')])

And because g was skipped in the 5 modification file, I know the issue is not e itself, just the number of mods per entry.

Do you know why this happens?
How could more than 5 modifications per base be processed with pysam?

Of note, we work with RNA modifications, there are over 170 of biologically concurring RNA modifications, so this will be of interest to wider audience soon.
Already at this time there are models to call 8 modifications in RNA.

[jmarshall]

Thanks for the clear bug report. Do you have test SAM files you could attach? I could generate some containing the MM/ML fields you show, but if you already have such files that'd save some time.

[lpryszcz]

I couldn't attached .bam here, so uploaded them to https://public-docs.crg.es/enovoa/public/lpryszcz/bugs/pysam/1348/

btw, I have found a workaround, simply processing MM/ML tags corresponding to up to 5 mods at the time works fine

step = 5
...
        if nmods>5:
            modified_bases = {}
            if not a.tags: continue
            tags = {k: v for k, v in a.tags}
            if 'MM' not in tags: continue
            parts = tags['MM'].split(";")[:-1]
            mlx = len(tags['ML'])//len(parts)
            for i in range(0, len(parts), step):
                a.set_tag('MM', ";".join(parts[i:i+step])+";")
                a.set_tag('ML', tags['ML'][i*mlx:(i+step)*mlx])
                # store
                for k, v in a.modified_bases.items():
                    modified_bases[k] = v

And I use modified_bases instead of a.modified_bases.

[jmarshall]

Thanks, I've downloaded those.

There is an array inside the modified_bases function that is hard-coded to 5 items. Just enlarging it will likely fix this — for a while! But it'll be better to see if we can figure out dynamically how large the array needs to be (maybe by counting semicolons in MM or something like that) or do something similar to your workaround.

[lpryszcz]

thanks @jmarshall, you mean here, right?

Unfortunately, I'm not fluent in C, so can't help much implementing it...

Maybe it could be smart to store a list of expected modifications in a BAM header? This way no guessing would be needed.
Additionally, having long names for mods encoded somewhere in BAM header would be useful @marcus1487

[jmarshall]

That is indeed the array in question.

Don't worry, I have a fix almost ready. There is no need for anything additional in any headers for fixing this: bam_mods_recorded() tells you the maximum required array size.

[jmarshall]

This fix depended on PR #1362, which took a while. Both have been applied and now pushed. So this fix will appear in the next release in due course.