Where the ring begins
Many viruses and bacteria have DNA closed into a ring: it has neither beginning nor end. The ring was read twice. The first reading starts somewhere, the second at a completely different place on the same ring, and on top of that the instrument made mistakes in some letters in the second reading. You have to work out where in the first reading the second one begins.
Formally: you are given strings s and t of the same length n over the letters A, C, G, T. For a shift d (0 ≤ d < n) the number of matches is the number of positions j with tj = s(j + d) mod n. Find the shift with the largest number of matches.
Input
Two lines: s and t, each of length n (1 ≤ n ≤ 2 · 10⁵), over the letters A, C, G, T.
Output
Two numbers: the shift d with the largest number of matches and that number itself. If there are several such shifts, print the smallest.
Examples
Example 1
ACGTTGCAAG CAAGACGTAG
6 9
Limits
- Time · Python
- 5000 ms
- Time · Go
- 2000 ms
- Memory
- 256 MB
- Judging
- 7 tests
Walkthrough
Open it step by step: each step gives away a little more. You can jump straight to the solution, but a problem solved with one hint stays with you longer.
1Hint 1
Comparing the strings at one shift is n operations, and there are n shifts: 4 · 10¹⁰ in all. For the phage φX174 of 5386 letters this is still possible, for the phage lambda of 48,502 it no longer is.
2Hint 2
Matches of letters are hard to add up, but numbers are easy. Take one letter, say A, and replace the string with zeros and ones: 1 where there is an A. What does the sum of products of these ones count at shift d?
3The idea
For a letter c let ai = 1 if si = c, and bj = 1 if tj = c. Then Σj bj · a(j + d) mod n is the number of positions where at shift d both strings show the letter c. The sum over the four letters is the number of matches we want. And such a sum for all d at once is a circular cross-correlation, which the fast Fourier transform computes in O(n log n): transform both arrays, multiply the first by the conjugate of the second, transform back.
4Walkthrough
In Python with numpy this is a few lines per letter: irfft(rfft(a) * conj(rfft(b)), n) gives a correlation of length exactly n, and no zero padding is needed: numpy can do transforms of any length, and a circular correlation is just what is needed here, since this is a ring.
Our own FFT in Go works with lengths that are powers of two, so the circle is unrolled: s is written twice in a row, t is reversed, and the ordinary convolution at index d + n − 1 gives the matches at shift d. The transform length is a power of two of at least 3n. To avoid three transforms per letter, both arrays of one letter go into one complex number: s into the real part, t into the imaginary part. The imaginary part of the square (x + iy)² is 2xy, which is exactly the product needed. The squares for the four letters are added up, and one inverse transform is done.
The second and third tests are real phages from the NCBI database: φX174, whose DNA is circular, and lambda, whose DNA closes into a ring in the cell by its sticky ends. Even with 30% errors in the fourth test the correct shift stands out: random matches give about a quarter of the positions, and the correct shift 70%.
4+Solution
5Complexity
Five transforms of length 2²⁰ in Go or twelve transforms of length n in Python: O(n log n). The correlation values are at most 2 · 10⁵, and rounding after the inverse transform recovers them exactly.
6Common mistakes
- All shifts head-on. Passes φX174 and fails on lambda: 48,502² ≈ 2.4 · 10⁹ comparisons.
- Finding the start by the first letters. Searching s for the first twenty letters of t is fast, but if they contain an instrument error there is no exact match at all, and with repeats in the genome there are several matches.
- An ordinary convolution instead of a circular one. Without doubling s, the tail of the second reading that runs past the end of the first is compared with nothing.
- Ties. The string ACGTACGT… matches itself at shifts 0, 4, 8, …, and the smallest is needed; the fifth test is about this.