ambiguous calls in nanopolish

Question

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?

Answer 1

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].num_reads += 1
    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.