Coordinate numbering between forward (+) and reverse (-) strand in Tombo

Question

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'

Answer 1

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.