Where the reads land
A reference genome is stored as its Burrows–Wheeler transform, as in the previous problem: a character $, smaller than all letters, is appended to a text T of the letters A, C, G, T, the cyclic rotations are sorted, and their last letters are written out, giving the string L. A sequencer has read the DNA in short pieces, reads. All reads are taken from the same strand as T, but some contain errors.
For each read, find where it lands. If the read occurs in T letter for letter, print all such places. If not, print all places where it matches T in all letters but one. If there are none of those either, the read does not land.
The tests include the genomes of phages λ, φX174 and three related phages in a row, human mitochondria, random letters and many copies of one piece. In the λ test the reads are taken from the laboratory strain with its four substitutions, as in the article.
Input
The first line is L: the letters A, C, G, T and one character $, of length 2 to 2·105 + 1. The second is the number of reads q (1 ≤ q ≤ 3000). Then q lines, reads of the letters A, C, G, T of length 1 to 40.
Output
For each read, a line: the number of places k followed by the places themselves in increasing order. A place is the number of the read's first letter in T, letters numbered from one. If the read does not land, the line is 0.
Examples
Example 1
TGTA$AAACCCGG 4 ACG CGA TTT GTACGT
3 1 5 9 1 6 0 2 3 7
Example 2
CTCGGTA$ATTAA 4 GATTAC CCCC GATTCC A
1 1 0 1 1 4 2 5 7 9
Limits
- Time · Python
- 8000 ms
- Time · Go
- 1500 ms
- Memory
- 256 MB
- Judging
- 10 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
Trying all places in the text for each read is 2·105 · 3000 · 40 comparisons, far too many. Exact occurrences come from backward search. How can it be allowed one error?
2Hint 2
Backward search reads a read from the end and picks a letter at each step. What if at one of the steps you take not the read's letter but another one, and remember that the error has been spent? How many such paths stay alive if a path is dropped as soon as its range is empty?
3The idea
First an ordinary backward search; if the range is non-empty, its rows are the exact occurrences. Otherwise backtracking: the state is (how many letters of the read remain to be read, the range, whether the error has been spent). From a state you can go by the read's letter, and if the error is not spent, also by the three other letters, spending it. An empty range ends the branch. The ranges reached after reading the whole read are rows where the read stands with one substitution; there is exactly one substitution, because there are no exact occurrences.
It remains to turn table rows into places in the text. The text itself is not needed for this: LF steps from row 0, which holds the rotation beginning with $, go through the rows of positions n − 1, n − 2, …, 0 one at a time, and the whole suffix array is written down along the way.
4Walkthrough
Keep rank counters for the five characters for all rows, so each step costs O(1). The suffix array is n numbers filled in by n LF steps. The backtracking is conveniently written with a stack rather than recursion. The places from all ranges found are sorted.
In the first sample the text is ACGTACGAACGT: ACG stands exactly at positions 1, 5 and 9; CGA exactly at position 6, and places with one substitution, like CGT at position 2, are no longer printed. TTT does not land. GTACGT does not occur exactly, and with one substitution it lands at positions 3 (GTACGA) and 7 (GAACGT).
4+Solution
5Complexity
Preparation is O(n). A read of length m costs O(σ·m2) steps in the worst case: for each place of the substitution, three letters and the rest of the read. In practice a branch with a substitution dies within a few steps: in a genome of 2·105 letters almost any ten random letters in a row do not occur. The output is as many places as were found.
6Common mistakes
- Printing both exact places and places with a substitution. If a read landed exactly, places with one substitution are not needed; the first sample catches this.
- Spending the error on the same letter. The branch “the read's letter, but with the error spent” gives the same rows as without spending, and places get doubled.
- Recursion in Python. The depth is only 40, but a function call in Python is expensive; a stack of tuples is faster.
- A pass over the text. Restoring the text and checking every place for each read is correct, but already the λ test does not fit in the time.