Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-01-11T23:34:52.743Z
Answers.Text: 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?
Answers.UserDisplayName: aerijman
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-01-11T19:44:10.537Z
Answers.LastActivityDate: 2021-01-11T19:44:10.537Z
Answers.UserDisplayName: Maximilian Press
Body: <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>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-01-11T15:40:52.427Z
Comments.Text: What's the question here? If there is not enough support for a position, it gets filtered out...
Comments.UserDisplayName: Wouter De Coster
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-01-11T19:04:50.050Z
Comments.Text: 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.
Comments.UserDisplayName: Maximilian Press
Tags: <nanopore><methylation>
Title: ambiguous calls in nanopolish
ThreadText: # ambiguous calls in nanopolish
  
  (Asked by: username_0 on 2021-01-10)
  
  I am converting scores into frequencies (from nanopolish results) using [Nanopore scripts](https://nanopolish.readthedocs.io/en/latest/quickstart_call_methylation.html).
  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 username_1 de username_1 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 username_1 de username_1'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](https://github.com/jts/nanopolish/blob/master/scripts/calculate_methylation_frequency.py):
  
  ```
  
  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.
  
  > **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). @username_1 de username_1, is that correct?  
  
  ---
  
  