How to output detailed information about the methylation sites occurring in each bin?

[luda99]

What should be done to enable nanomitf to output information about each methylation site, methylation type, and methylation score? For now, the file "bin-motifs.tsv" only contains the quantity. I think since we can output the quantity, it must be possible to output detailed methylation site information for each Bin as well.

[SebastianDall]

Hi @luda99

Do you mean a pileup type file where each motif is mapped where you can see the n_modified and v_valid_cov?

If so you can run the epimetheus cli tool (https://github.com/SebastianDall/epimetheus). You can also run this via the python bindings, which should already be installed via nanomotif conda env.

Python:

import os
from epymetheus import epymetheus
from epymetheus.epymetheus import MethylationOutput


def methylation_pattern_at_positions():
    pileup = os.path.join("path/to/", "geobacillus-plasmids.pileup.bed")
    assembly = os.path.join("path/to/", "geobacillus-plasmids.assembly.fasta")

    outfile = os.path.join("path/to/","out.tsv")

    epymetheus.methylation_pattern(
        pileup,
        assembly,
        motifs = ["GATC_a_1", "GATC_m_3", "RGATCY_a_2"], # Motifs from bin-motfis.tsv
        output = str(outfile),
        threads = 1,
        output_type=MethylationOutput.Raw
    )

rust cli:

epimetheus methylation-pattern contig -p /path/to/pileup.bed -a /path/to/assembly.fasta -m GATC_a_1 RGATCY_a_2 ... --threads 1 --output-type raw --output /path/to/out.tsv

The output will look like so:

contig	start	strand	motif	mod_type	mod_position	n_modified	n_valid_cov
contig_2	281	+	GATC	a	1	7	8
contig_2	282	-	GATC	a	1	1	3
contig_2	513	+	GATC	a	1	8	9
contig_2	514	-	GATC	a	1	3	4
contig_2	621	+	GATC	a	1	8	9
contig_2	622	-	GATC	a	1	5	5
contig_2	714	+	GATC	a	1	10	10
contig_2	715	-	GATC	a	1	5	6
contig_2	756	+	GATC	a	1	11	11

Where you can see n_modified and n_valid cov of each GATC_a_1 observation in the contig denoted by "start". In essence this is a pileup filter.

I hope this helps.

Kind regards,
Sebastian

[luda99]

Thank you for your prompt reply. In other words, I would like to know exactly which base positions are represented by the 'n_mod' and 'n_nomod' values output in the bin-moitf.tsv file?
For example: reference motif mod_position mod_type n_mod n_nomod motif_type motif_complement mod_position_complement n_mod_complement n_nomod_complement
A12C_bin.10.fa CNTCCGGANA 4 21839 218 41 non-palindrome
A12C_bin.10.fa contains only 259 CNTCCGGANA motifs? Where are the 218 mods exactly? Where are the 41 non-motifs exactly?
Because later on, I plan to compare the differentially methylated positions of the same MAG in different sample metagenomes. Of course, I have confirmed using instrain that this MAG is the same in different samples.
Thank you again. @SebastianDall

[SebastianDall]

Ah, thank you for the clarification.

I would:

1. extract the contig ids of A12C_bin.10.fa
2. run:
   pileup = os.path.join("path/to/", "pileup.bed.gz")
   assembly = os.path.join("path/to/", "assembly.fasta")

This will be much faster if the pileup is bed.gz! If not it is okay but may take a long time depending on the pileup size.

df = epymetheus.methylation_pattern(
pileup,
assembly,
contigs = ["contigs_in_A12C_bin.10.fa"],
motifs = ["CNTCCGGANA_21839_4"], # Motifs from bin-motfis.tsv
# output = str(outfile), # optional specify a output file.
threads = 1,
output_type=MethylationOutput.Raw
)

3. The above returns a polars dataframe.
   To find which are deemed methylated just find the positions which are >=70% methylated.

df = df.with_columns(
((col("n_modified") / col("n_valid_cov")) >= 0.70).alias("is_modified")
)

This should give you the positions nanomotif deems as methylated. The 41 positions that you mention could be due to strain difference where the fraction_modified just dips under 0.7 rather than it being 0 and 100 % methylated. So I would check the fraction modified at these positions.

There may be a small discrepancy between how Nanomotif calculates n_valid_cov and epimetheus but it should hopefully be negligable if there is one.

Hope that helps.

[luda99]

Thank you. I will try your method. @SebastianDall