Answers.Body: <p>When <strong>Tombo re-squiggle</strong> and aligns a read it adds two attributes to it: <strong>mapped_start</strong> and <strong>mapped_end</strong>. Tombo also gives chromosomes and strands.</p>
  <p>Imagine you are looking for <strong>some value of a signal[s] in position x to y on a chromosome on both strands.</strong> In this range, you probably have reads that span the positive strand (forward) and the reads that cover the negative strand (reverse).</p>
  <p>Now let a given read that has mapped to a positive strand and covers the entire [x, y], have attributes: <strong>pos_mapped_start</strong> and <strong>pos_mapped_end</strong>. Similarly let the reads mapped to negative strands covering [x, y] have attributes: <strong>neg_mapped_start</strong> and <strong>neg_mapped_end</strong>.</p>
  <p>For reads mapped to + strand:</p>
  <pre><code>start_idx = x- pos_mapped_start
  end_idx = y - pos_mapped_start + 1
  required_signal= read_signal[start_idx: end_idx]
  </code></pre>
  <p>For reads mapped to - strand:</p>
  <pre><code>start_idx = neg_mapped_end - y- 1
  end_idx = neg_mapped_end - x
  required_signal= read_signal[start_idx: end_idx]
  </code></pre>
  <ul>
  <li>Note that x and y are <strong>inclusive</strong> and the length of the required signal = y - x + 1.</li>
  </ul>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-09-17T13:43:01.283Z
Answers.LastActivityDate: 2022-09-17T13:43:01.283Z
Answers.UserDisplayName: Marjan
Body: <p>My question is related to the output of Tombo re-squiggling algorithm on nanopore data.
  Given the position/coordinate on one strand (+ or -), I want to find the coordinate of its exact location on the opposite strand.
  For example, let's say I want to find data/signal on specific coordinates on both strands on Tombo path after re-squiggling. Let's say position 13936 - 14036 on genome on + and - negative strands. I wrote the following code but the sequence obtained from Tombo is not similar at all. I want to know how I can find the coordinate corresponding to a given position on the opposite strand.</p>
  <pre><code># one read from the + strand
  mapped_start, mapped_end, read_id, fast5_path, _, _ = reads_per_win_pos.loc[rd_id, :]
  &gt;&gt;&gt; start
  13936
  &gt;&gt;&gt; end
  13942
  &gt;&gt;&gt; mapped_start, mapped_end, read_id, fast5_path, _, _ = reads_per_win_pos.loc[rd_id, :]
  &gt;&gt;&gt; mapped_start
  11205
  &gt;&gt;&gt; mapped_end
  20084
  &gt;&gt;&gt; f = h5py.File(fast5_path, 'r')
  &gt;&gt;&gt; tombo_events = list(f['Analyses/RawGenomeCorrected_000/BaseCalled_template/Events'])
  &gt;&gt;&gt; forward_seq = ''.join([e[4].decode() for e in tombo_events])
  &gt;&gt;&gt; start_idx = start - mapped_start
  &gt;&gt;&gt; start_idx
  2731
  &gt;&gt;&gt; win_seq = seq[start_idx: start_idx + 100]
  &gt;&gt;&gt; forward_win_seq = forward_seq[start_idx: start_idx + 100]
  &gt;&gt;&gt; # one other overlapping read on the other strand (-)
  ...
  &gt;&gt;&gt; mapped_start, mapped_end, read_id, fast5_path, _, _ = reads_per_win_neg.loc[rd_id, :]
  &gt;&gt;&gt; mapped_start
  10979
  &gt;&gt;&gt; mapped_end
  20091
  &gt;&gt;&gt; f = h5py.File(fast5_path, 'r')
  &gt;&gt;&gt; tombo_events = list(f['Analyses/RawGenomeCorrected_000/BaseCalled_template/Events'])
  &gt;&gt;&gt; reverse_seq = ''.join([e[4].decode() for e in tombo_events])
  &gt;&gt;&gt; start_idx = start - mapped_start
  &gt;&gt;&gt; start_idx
  2957
  &gt;&gt;&gt; reverse_win_seq = reverse_seq[start_idx: start_idx + 100]
  </code></pre>
  <p>Now the forward_win_seq and reverse (or even the reverse complement) of it are not similar:</p>
  <pre><code>&gt;&gt;&gt; forward_win_seq
  'GAGACGGAGCAGACCCATCTGCTACTGCCCTTTCTATAATAACTAAAGTTAGCTGCCCTGGACTATTCACCCCCTAGTCTCAATTTAAAAAGATCCCCAT'
  &gt;&gt;&gt; reverse_win_seq
  'GGGCCAGGTGCCTGAGATCCATGTTCCATCCTACCTGCCTGACCTGCCCGGCATTGCCAACGACCTCATGTACAGTGCCGACCTGGGCCCCGGCATTGCC'
  &gt;&gt;&gt; complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}
  &gt;&gt;&gt; com_seq = ''.join([complement[base] if base in complement.keys() else base for base in reverse_win_seq])
  &gt;&gt;&gt; com_seq[::-1]
  'GGCAATGCCGGGGCCCAGGTCGGCACTGTACATGAGGTCGTTGGCAATGCCGGGCAGGTCAGGCAGGTAGGATGGAACATGGATCTCAGGCACCTGGCCC'
  </code></pre>
  
Tags: <nanopore><read-mapping>
Title: Coordinate numbering between forward (+) and reverse (-) strand in Tombo
ThreadText: # Coordinate numbering between forward (+) and reverse (-) strand in Tombo
  
  (Asked by: username_0 on 2022-09-16)
  
  My question is related to the output of Tombo re-squiggling algorithm on nanopore data.
  Given the position/coordinate on one strand (+ or -), I want to find the coordinate of its exact location on the opposite strand.
  For example, let's say I want to find data/signal on specific coordinates on both strands on Tombo path after re-squiggling. Let's say position 13936 - 14036 on genome on + and - negative strands. I wrote the following code but the sequence obtained from Tombo is not similar at all. I want to know how I can find the coordinate corresponding to a given position on the opposite strand.
  
  ```
  # one read from the + strand
  mapped_start, mapped_end, read_id, fast5_path, _, _ = reads_per_win_pos.loc[rd_id, :]
  >>> start
  13936
  >>> end
  13942
  >>> mapped_start, mapped_end, read_id, fast5_path, _, _ = reads_per_win_pos.loc[rd_id, :]
  >>> mapped_start
  11205
  >>> mapped_end
  20084
  >>> f = h5py.File(fast5_path, 'r')
  >>> tombo_events = list(f['Analyses/RawGenomeCorrected_000/BaseCalled_template/Events'])
  >>> forward_seq = ''.join([e[4].decode() for e in tombo_events])
  >>> start_idx = start - mapped_start
  >>> start_idx
  2731
  >>> win_seq = seq[start_idx: start_idx + 100]
  >>> forward_win_seq = forward_seq[start_idx: start_idx + 100]
  >>> # one other overlapping read on the other strand (-)
  ...
  >>> mapped_start, mapped_end, read_id, fast5_path, _, _ = reads_per_win_neg.loc[rd_id, :]
  >>> mapped_start
  10979
  >>> mapped_end
  20091
  >>> f = h5py.File(fast5_path, 'r')
  >>> tombo_events = list(f['Analyses/RawGenomeCorrected_000/BaseCalled_template/Events'])
  >>> reverse_seq = ''.join([e[4].decode() for e in tombo_events])
  >>> start_idx = start - mapped_start
  >>> start_idx
  2957
  >>> reverse_win_seq = reverse_seq[start_idx: start_idx + 100]
  
  ```
  
  Now the forward\_win\_seq and reverse (or even the reverse complement) of it are not similar:
  
  ```
  >>> forward_win_seq
  'GAGACGGAGCAGACCCATCTGCTACTGCCCTTTCTATAATAACTAAAGTTAGCTGCCCTGGACTATTCACCCCCTAGTCTCAATTTAAAAAGATCCCCAT'
  >>> reverse_win_seq
  'GGGCCAGGTGCCTGAGATCCATGTTCCATCCTACCTGCCTGACCTGCCCGGCATTGCCAACGACCTCATGTACAGTGCCGACCTGGGCCCCGGCATTGCC'
  >>> complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}
  >>> com_seq = ''.join([complement[base] if base in complement.keys() else base for base in reverse_win_seq])
  >>> com_seq[::-1]
  'GGCAATGCCGGGGCCCAGGTCGGCACTGTACATGAGGTCGTTGGCAATGCCGGGCAGGTCAGGCAGGTAGGATGGAACATGGATCTCAGGCACCTGGCCC'
  
  ```
  
  **Tags:** nanopore, read-mapping
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_0 on 2022-09-17. Score: 2 | Accepted answer
  
  When **Tombo re-squiggle** and aligns a read it adds two attributes to it: **mapped\_start** and **mapped\_end**. Tombo also gives chromosomes and strands.
  
  Imagine you are looking for **some value of a signal[s] in position x to y on a chromosome on both strands.** In this range, you probably have reads that span the positive strand (forward) and the reads that cover the negative strand (reverse).
  
  Now let a given read that has mapped to a positive strand and covers the entire [x, y], have attributes: **pos\_mapped\_start** and **pos\_mapped\_end**. Similarly let the reads mapped to negative strands covering [x, y] have attributes: **neg\_mapped\_start** and **neg\_mapped\_end**.
  
  For reads mapped to + strand:
  
  ```
  start_idx = x- pos_mapped_start
  end_idx = y - pos_mapped_start + 1
  required_signal= read_signal[start_idx: end_idx]
  
  ```
  
  For reads mapped to - strand:
  
  ```
  start_idx = neg_mapped_end - y- 1
  end_idx = neg_mapped_end - x
  required_signal= read_signal[start_idx: end_idx]
  
  ```
  
  * Note that x and y are **inclusive** and the length of the required signal = y - x + 1.
  
  
  ---
  
  