Using your own modkit-generated pileups with plot_enrichment_profile.by_modification() instead of BAM input to parse_bam.pileup()

[K1999ban]

Hi,
I generated pileup files using modkit, for both 5mC (CG,0) and 6mA (A,0), from BAM files that contain MM and ML tags.

I then reformatted the modkit-generated pileups to match the output of dimelo_v2.parse_bam.pileup() — including removing extra columns, bgzip compression, and tabix indexing.

However, when I try to use them as input to plot_enrichment_profile.by_modification(), like this:

plot_enrichment_profile.by_modification(
mod_file_name='modkit_pileup.tsv.gz',
regions=CTCF_BED,
window_size=1000,
motifs=['CG,0', 'A,0'],
smooth_window=50
)
I get a ValueError: Modification positions out of expected window range

Can the plotting functions in dimelo_v2 work with pileups generated by different modkit versions, as long as the format matches? Or is there something specific in DiMeLo's own pileup output that the plot function needs? Because I run modkit_v0.2.4 for Dimelo V2 when I generate the pileups from BAM file normally.

Also, is there a way to filter by coverage within DiMeLo v2 when generating pileups?

For example, when using modkit, the pileup output (BEDmethyl .bed) has 18 columns, and column 10 corresponds to n_valid_cov. I'd be interested in filtering for positions where n_valid_cov >= 5 to retain only high-confidence base calls.

Is there currently a way to apply such a filter in DiMeLo v2 before or during plotting?
Or should this filtering be done externally on the pileup file before passing it to plot_enrichment_profile.by_modification()?

Thanks a lot!

[OberonDixon]

Wow, sorry we missed this. We were focused for a while on new assay development work and merging the dimelo_v2 code back into the main dimelo repo: https://github.com/streetslab/dimelo

If you're still struggling with anything or have feature requests you can post there. Currently we do not support filtering on loading from the pileup bedmethyl files (or equivalently formatted tsvs) directly, but if you load the pileup vectors using the load_processed.pileup_vectors_from_bedmethyl function then the numerator (modified bases) and denominator (number of valid reads for the site) are both available so filtering could be done at that step. It would not be difficult to integrate into the main plotters we would just need to understand the scope of the use case to design something reasonable.

[K1999ban]

Thanks a lot for the detailed response! That clarifies it.

I’ll try loading the modkit-generated pileups through
load_processed.pileup_vectors_from_bedmethyl()
and perform the coverage filtering at that step using the denominator values.

Just to confirm- I will filter the pileup vectors before passing them to plot_enrichment_profile.by_modification(), that could avoid the “Modification positions out of expected window range”.

Also, it would be great if future versions of dimelo_v2 supported an optional coverage cutoff argument (e.g. min_cov) within the plot functions- it would simplify high-confidence analyses.

Thanks again for the clarification!

[OberonDixon]

That's an interesting feature request, thanks! By the way, I think I actually didn't directly answer the modkit version question: we plan to very soon add support for the newest modkit versions, but currently they are not supported because the bedmethyl format was changed and some stuff on our end is hardcoded to the old format. When the change is implemented, it'll be on the main dimelo repo https://github.com/streetslab/dimelo, as this one is no longer being updated.

In terms of implementing your own filtering: if want to do this filtering across many regions, counting only coordinates with sufficient depth in each specific region, then you should actually use load_processed.regions_to_list with function_handle=pileup_vectors_from_bedmethyl and the appropriate kwargs. This will return a tuple of vectors for each individual region individually, which you can then filter.

def regions_to_list(
    function_handle,
    regions,
    window_size: int | None = None,
    quiet: bool = True,
    cores: int | None = None,
    split_large_regions: bool = False,
    **kwargs,
):
    """
    User-facing function.

    Run any standard load_processed pileup or extract loader loading each region from the region
    specifier into a new element of a list.

    Args:
        function_handle: the loader function you want to run.
        regions: the region specifier. Typically we expect to get many regions for this function, in the form of a list
            of strings or bed file paths. regions_to_list will run across all of these one-by-one returning a separate
            function return for each independent region.
        window_size: window around centers of regions, defaults to None
        quiet: disables progress bars
        cores: CPU cores across which to parallelize processing. Default to None, which means all available.
        split_large_regions: if True, regions will be run sequentially in parallelized chunks. If False,
            each individual region's chunks will be run sequentially but there will be parallelization across

However, if you only care about whether there are enough counts total in terms of coverage for the given positions, then the direct pileup_vectors_from_bedmethyl call with the set of regions can give you that all summed up so you can filter.

To pass custom data (e.g. filtered in this way) to the enrichment plotter, you actually want to send a list of trace_vector values to plot_enrichment_profile.make_enrichment_profile_plot, with corresponding sample_names.

def make_enrichment_profile_plot(
    trace_vectors: list[np.ndarray],
    sample_names: list[str],
    **kwargs,
) -> Axes:
    """
    Plot the given enrichment profile traces.

    This helper function can be useful during plot prototyping, when repeatedly building plots from the same data.
    The first argument should be the output of get_enrichment_profiles().

    Args:
        trace_vectors: list of enrichment profile traces
        sample_names: list of names to use for labeling traces in the output; legend entries
        kwargs: other keyword parameters passed through to utils.line_plot

    Returns:
        Axes object containing the plot
    """

You'd want to calculate the mod/valid fraction using code something like this:

                nans_everywhere = np.full_like(
                    modified_base_counts, np.nan, dtype=float
                )
                trace = np.divide(
                    modified_base_counts,
                    valid_base_counts,
                    out=nans_everywhere,
                    where=valid_base_counts != 0,
                )