0
$\begingroup$

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?

$\endgroup$
2
  • 1
    $\begingroup$ What's the question here? If there is not enough support for a position, it gets filtered out... $\endgroup$ Jan 11, 2021 at 15:40
  • $\begingroup$ 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. $\endgroup$ Jan 11, 2021 at 19:04

1 Answer 1

1
$\begingroup$

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.
$\endgroup$
1
  • $\begingroup$ 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? $\endgroup$
    – aerijman
    Jan 11, 2021 at 23:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.