<p>Assuming that it is the postprocessing script <code>calculate_methylation_frequency.py</code> 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 <code>LLR &gt; T * N</code>, but also the statistical threshold <code>LLR &gt; 0</code> (see below).</p>
<p>Specifically, here is the code in the <a href="https://github.com/jts/nanopolish/blob/master/scripts/calculate_methylation_frequency.py" rel="nofollow noreferrer">script</a>:</p>
<pre><code>
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 &gt; 0:
        sites[key].called_sites_methylated += num_called_cpg_sites

</code></pre>
<p>That is the function that decides whether to add a site to output, which depends on <code>is_methylated</code>, which is set here:</p>
<pre><code>    # [skip several lines of code]
        num_sites = int(record['num_motifs'])
        llr = float(record['log_lik_ratio'])

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

        # [MP NOTE: SKIP IF LLR &lt;= 0]
        is_methylated = llr &gt; 0

        # [skip several lines of code]

        update_call_stats(key, num_sites, is_methylated, sequence) 
</code></pre>
<p>It's a little messy to read because the <code>sites</code> dictionary is a global outside of <code>update_call_stats()</code>.</p>
<p>So overall, it seems like <strong>if you want sites included that are below this threshold, you will have to either:</strong></p>
<ol>
<li>use an earlier stage of processing (output of <code>nanopolish call-methylation</code> I believe according to the docs) or</li>
<li>alter the script to change the heuristic thresholds.</li>
</ol>

CC BY-SA 4.0
2021-01-11T23:34:52.743Z
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?
aerijman
CC BY-SA 4.0
2021-01-11T19:44:10.537Z
2021-01-11T19:44:10.537Z
Maximilian Press
<p>I am converting scores into frequencies (from nanopolish results) using <a href="https://nanopolish.readthedocs.io/en/latest/quickstart_call_methylation.html" rel="nofollow noreferrer">Nanopore scripts</a>.
The script filters out lots of rows on the premise that reads where <code>LLR &lt; T * N</code> are ambiguous.</p>
<p>where:</p>
<pre><code>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.
</code></pre>
<p>What? Why? Where?</p>

CC BY-SA 4.0
2021-01-11T15:40:52.427Z
What's the question here? If there is not enough support for a position, it gets filtered out...
Wouter De Coster
CC BY-SA 4.0
2021-01-11T19:04:50.050Z
I'm assuming that you are using the workflow defined here: https://nanopolish.readthedocs.io/en/latest/quickstart_call_methylation.html . 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.
Maximilian Press
<nanopore><methylation>
ambiguous calls in nanopolish
# ambiguous calls in nanopolish

(Asked by: username_0 on 2021-01-10)

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?

**Tags:** nanopore, methylation
> **Comments:**
> **username_1**: What's the question here? If there is not enough support for a position, it gets filtered out...  
> **username_2**: I'm assuming that you are using the workflow defined here: https://nanopolish.readthedocs.io/en/latest/quickstart_call_methylation.html . 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.  

---

## 1 answers

### Answer by username_2 on 2021-01-11. Score: 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:

use an earlier stage of processing (output of nanopolish call-methylation I believe according to the docs) or
alter the script to change the heuristic thresholds.

> **Comments:**
> **username_0**: 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?  

---

