2
$\begingroup$

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'
$\endgroup$

1 Answer 1

2
$\begingroup$

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.
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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