request enhancements to modBAM parsing

[sathish-t]

Hello Dr. Vollger,

I am writing a rust package to deal with modification data in a different context than fiberseq. Since fibertools parses Mm and Ml tags in modBAM data and is open source, I am thinking of importing the function my_mm_ml_parser and the like from src/utils/basemods.rs in my package.

This function doesn't cover all Mm/Ml tag formats, so I was wondering if you have plans to make it more general; please see the two examples below. I mean I could take a crack at it with some guidance and do a pull request here, or just take your code with an acknowledgement and put it in my package and try to modify it there. My other options are (1) to use rust_htslib which doesn't do an out of the box reference to query mapping of modification data, which you do as far as I can tell, or (2) to use modkit from ONT which is not under an MIT license.

(1) format A: MM:Z:C+mh,5,12; ML:B:C,204,26,89,130 is equivalent to format B: ‘MM:Z:C+m,5,12;C+h,5,12; ML:B:C,204,89,26,130’ As far as I can tell, you support the B format but not the A. I've not seen any example of software using the A format, but maybe in the future some technology calling two types of methylation on the C will use this.

(2) ‘MM:Z:C+m,5,12;C+h,5,12; ML:B:C,204,89,26,130’, this tag technically means any cytosine on the molecule not mentioned in the tag is to be regarded as unmodified. But your code treats any cytosine not mentioned in the tag as "missing" rather than unmodified if I am not mistaken.

Regards
Sathish

[mrvollger]

You are right on point 1 for sure. I was not thinking about 2 character codes. Any chance you are interested in writing a PR?

If not you can wait on me (might be a while) or use rust htslib + my bamlift module. (I think)

Point 2 I am not sure I understand. I should add code to deal with the ?/. format. But more directly to your question, and in the "." case, what change are you suggesting? I am not sure how else I'd represent them.

Best,
Mitchell

[sathish-t]

Hi Mitchell,

I'll try writing a PR. I am quite new to Rust so it will take me some time!

Point 2, I think I misunderstood the code.. I thought my_mm_ml_parser outputs both modified and unmodified positions along a sequence but I now think it outputs only the modified positions. Is this correct? Then I don't really have an issue here.

Regards
Sathish

[mrvollger]

Yep that code chunk only outputs modified positions.

[sathish-t]

Hi Mitchell,

Would you be interested in a new option to strip-basemods that would remove modifications that fall outside a given bed file? Like ft strip-basemods --include-bed $bed_file $input_bam $output_bam where $bed_file can be one of the following two formats:

chr1\t100\t200
chr2\t300\t400

in which case modifications outside these reference positions and on unmapped reads will get removed, or a format

chr1\t100\t200\tsome_read_name
chr2\t300\t400\tsome_other_read_name

where only modifications on some_read_name and some_other_read_name within these intervals are retained and everything else is dropped?

I feel like if I want to calculate an auto correlation function in fiberseq but only for a set of genomic regions, or only for specific reads, this function might come in handy. I have an independent requirement for this from some datasets unrelated to fiber-seq but which are in the mod BAM format. Samtools can filter by read or by region but it cannot filter modification data using these criteria.

I can edit the iterator in subcommands/strip_basemods.rs and any other file as needed to add these functions when I have the time and open a pull request if you are interested.

Regards
Sathish

[mrvollger]

I would be interested! Sounds really useful.

In terms of where you add the code, you might also consider looking at our filter expression module that has access to filtering data in any ft subcommand.

It would be a little harder than what you suggest but would work anywhere in ft.

Alternatively if you make the PR in strip-base mods I could extend it to be a filter expression.

Thanks!
Mitchell

[mrvollger]

I would be interested! Sounds really useful.

In terms of where you add the code, you might also consider looking at our filter expression module that has access to filtering data in any ft subcommand.

It would be a little harder than what you suggest but would work anywhere in ft.

Alternatively if you make the PR in strip-base mods I could extend it to be a filter expression.

Thanks!
Mitchell