# ambiguous calls in nanopolish

I am converting scores into frequencies (from nanopolish results) using Nanopore scripts. The script filters out lots of rows on the premise that reads where LLR < T * N are ambiguous.

where:

LLR = log-likelihood ratio (log-likelihood methylated / log-likelihood unmethylated)
T   = arbitrary threshold. Above that value the C is methylated.
N   = number of Cs (CpG) in the fragment analyzed.


What? Why? Where?

• What's the question here? If there is not enough support for a position, it gets filtered out... Jan 11, 2021 at 15:40
• I'm assuming that you are using the workflow defined here: nanopolish.readthedocs.io/en/latest/… . Can you give more information about where you are in the process of the "Nanopore scripts?" I am guessing that you are at calculate_methylation_frequency.py based on available information, and in that case it is as simple as Wouter de Coster says. Jan 11, 2021 at 19:04

Assuming that it is the postprocessing script calculate_methylation_frequency.py that you are using, user Wouter de Coster's comment is correct that the methylation event is filtered out at the frequency estimation stage if it has a low LLR. It will have a low LLR if it is considered not confident enough to report. This means not only LLR > T * N, but also the statistical threshold LLR > 0 (see below).

Specifically, here is the code in the script:


def update_call_stats(key, num_called_cpg_sites, is_methylated, sequence):
if key not in sites:
sites[key] = SiteStats(num_called_cpg_sites, sequence)

sites[key].called_sites += num_called_cpg_sites
if is_methylated > 0:
sites[key].called_sites_methylated += num_called_cpg_sites



That is the function that decides whether to add a site to output, which depends on is_methylated, which is set here:

    # [skip several lines of code]
num_sites = int(record['num_motifs'])
llr = float(record['log_lik_ratio'])

# Skip ambiguous call
if abs(llr) < args.call_threshold * num_sites:
# [MP NOTE: SKIP if LLR < T * N]
continue
sequence = record['sequence']

# [MP NOTE: SKIP IF LLR <= 0]
is_methylated = llr > 0

# [skip several lines of code]

update_call_stats(key, num_sites, is_methylated, sequence)


It's a little messy to read because the sites dictionary is a global outside of update_call_stats().

So overall, it seems like if you want sites included that are below this threshold, you will have to either:

1. use an earlier stage of processing (output of nanopolish call-methylation I believe according to the docs) or
2. alter the script to change the heuristic thresholds.
• Thanks! It wasn't easy to follow the script. I understood that the whole fragment is discarded if it can be defined which CpG is methylated (if at all). @Wouter de Coster, is that correct? Jan 11, 2021 at 23:34