Преобразование Барроуза — Уилера

ДНК фага λ, которую продают лабораториям, режут ферментом HindIII и получают восемь кусков известной длины — по ним сверяют длины других кусков ДНК в разделении на геле. Но если искать слово AAGCTT, где режет HindIII, в эталонном геноме λ из базы NCBI, мест окажется шесть, а кусков — семь. Лабораторный λ — особый штамм, cI857 ind1 Sam7, и от эталона он отличается четырьмя буквами из 48 502, одна из которых и создаёт седьмое место. Как найти эти буквы? Так же, как ищут отличия любого генома: разрезать ДНК на короткие куски, прочитать их — это риды — и для каждого найти место в эталоне. Ниже двадцать тысяч смоделированных ридов этого штамма, и все их ищет указатель, построенный на преобразовании Барроуза — Уилера.

Под лупой — позиция 37 589. В эталоне там C, а у ридов, которые легли на это место, — T. Буквы вокруг читаются как AAGCTT: это и есть седьмое место HindIII. Кнопка «к следующему отличию» обходит все места, где риды не согласны с эталоном. При ошибке ридов 1 % из двадцати тысяч 12 005 ложатся точно, 6165 — с одной заменой, 1830 не находятся, а отличий четыре, и все четыре — настоящие замены штамма. Поставьте ошибку 3 %. Каждый рид, который лёг на замену, уже потратил на неё свою единственную разрешённую ошибку, и рид с ещё одной ошибкой не ложится вовсе. Не находятся уже 9047 ридов, на позицию 37 589 ложатся только три, и сцена видит две замены из четырёх: поиск с одной ошибкой тянет картину к эталону.

Идея

Отсортируем все суффиксы текста, как в суффиксном массиве, и выпишем букву, стоящую перед каждым, — это преобразование Барроуза — Уилера. Одна эта строка — та же длина, что у текста, — позволяет найти любое слово: читать его с конца и после каждой буквы сужать отрезок строк таблицы, которые начинаются с уже прочитанного хвоста. Каждый шаг — два подсчёта «сколько раз такая буква встретилась в строке до этого места».

Само преобразование и то, почему оно собирает одинаковые буквы в серии, разобраны в главе курса «Конкурс упаковки»; там же — как bzip2 его сжимает и как текст восстанавливается обратно. Здесь — другая сторона: преобразование как указатель.

Допишем к тексту T длины n знак $, который меньше всех букв и больше нигде не встречается, и выпишем все циклические сдвиги T$ в алфавитном порядке. Получится таблица из n + 1 строк. Её первый столбец F — буквы текста по алфавиту, последний L — и есть преобразование. Так как $ стоит в конце и он самый маленький, порядок сдвигов совпадает с порядком суффиксов, и L[i] — буква перед суффиксом номер SA[i]: L[i] = T[SA[i] − 1], а для суффикса, который начинается с начала текста, — $. Для фага λ первые буквы L — GCGCGGAGAAAGGAG…: в 48 503 знаках преобразования 35 329 серий одинаковых букв. У ДНК одного генома серии короткие — повторов в ней мало.

Последний столбец ведёт в первый

Обозначим через C[c] число знаков T$, которые меньше c, считая и $: строки, начинающиеся с c, занимают в таблице отрезок с C[c] по C[c] + (число c) − 1. А через rank_c(i) — сколько раз буква c встречается в L[0..i), до строки i. Всё держится на одном наблюдении.

Шаг LF переводит строку, где текст начинается с позиции p, в строку, где он начинается с позиции p − 1: на букву левее. Если начать со строки 0 — там сдвиг «$T», и последняя буква — последняя буква текста, — и делать шаги LF, по столбцу L прочитаются все буквы текста справа налево. Это обратное преобразование за n шагов, и именно его надо найти самому в задаче главы курса.

Обратный поиск

Строки таблицы, которые начинаются с образца P, стоят подряд — они отсортированы. Найдём этот отрезок, читая P с конца.

Начинаем с отрезка [0, n + 1) — с ним начинаются все строки, — и прибавляем буквы образца по одной справа налево. После m шагов длина отрезка — число вхождений образца. Каждый шаг — два вызова rank, и если rank отвечает за постоянное время, образец считается за O(m) — независимо от длины текста. В листинге место EcoRI, GAATTC, ищется в геноме λ за 12 вызовов rank: по два на букву.

Остаётся сделать rank быстрым. Хранить для каждой буквы и каждой строки готовый ответ — это σ чисел на строку, больше самого текста. Поэтому ответы хранят через каждые 64 строки, а остаток досчитывают по L: не больше 63 букв. В настоящих реализациях блок L хранят битами, и досчёт — одна-две инструкции подсчёта единиц в машинном слове; для больших алфавитов есть вейвлет-деревья, где rank стоит O(log σ).

Здесь обратный поиск идёт по «Войне и миру». Слово «наташа» читается с конца. После «а» отрезок огромный: с этой буквы начинаются 183 707 суффиксов. После «ша» их 3908, после «аша» — 972, после «таша» — 827, и дальше число не меняется: всякое «таша» в романе — из имени «Наташа». Попробуйте одну букву «и» и сравните с длинной фразой: число вызовов rank зависит только от длины того, что вы ищете.

Главное в этой сцене — чего в ней нет. Браузер не загружал текст романа: только преобразование Барроуза — Уилера и номера позиций для каждой 256-й строки. Счёт, места на полосе и даже отрывки вокруг найденного он получил из этой одной строки. Файл с перестановкой и метками, сжатый gzip, весит 1,0 МБ, а сам роман, сжатый так же, — 1,4 МБ: серии одинаковых знаков сжимаются лучше текста. В перестановке романа 1 293 114 серий на 2 904 145 знаков, в среднем по 2,25 — у прозы серии длиннее, чем у ДНК фага.

Где стоит найденное

Отрезок говорит, сколько раз встречается образец, но не где. Место строки — это SA[i], а хранить весь суффиксный массив — n чисел по 4 байта, вчетверо больше текста. Поэтому хранят только каждую s-ю позицию: в листинге — позиции, кратные 32, вместе с номерами их строк. Чтобы найти позицию строки i, делают шаги LF, пока не попадут в отмеченную строку: каждый шаг сдвигает на букву левее, и если после k шагов попали в строку позиции q, то искомая позиция — q + k.

Так память и время меняются местами: при шаге 32 номера занимают в 32 раза меньше места, чем весь массив, а каждое место стоит до 31 шага. Сцена с романом отмечает каждую 256-ю строку — так файл меньше, — и там шагов в среднем около 256, но без гарантии; поэтому на полосе не больше трёх тысяч мест. Той же выборкой, только обратной — «строка для каждой 256-й позиции», — сцена читает отрывки: начинает с ближайшей отмеченной позиции справа от нужного куска и идёт LF налево, выписывая буквы из L.

Риды: поиск с одной ошибкой

Риды почти никогда не совпадают с эталоном буква в букву: машина ошибается, а у самого организма бывают свои замены — те, которые мы и ищем. Обратный поиск легко научить прощать замены. На шаге, где читается буква P[k], можно продолжить не только с ней, но и с каждой из трёх других букв, пометив, что единственная разрешённая замена потрачена. Получается перебор с возвратом по дереву вариантов, но дерево маленькое: ветка обрывается, как только отрезок пустеет, а в геноме из 48 тысяч букв почти любые десять случайных букв подряд не встречаются, так что ветки с заменой у конца рида гибнут через несколько шагов.

Листинг делает это с двумя тысячами ридов по 35 букв с обеих цепей генома, где каждая буква ошибочна с вероятностью 1 %. 1409 ридов легли точно, 499 — с одной заменой, 92 не нашлись: в них две ошибки и больше. Это и ожидалось: без ошибок остаётся 0,99³⁵ ≈ 70 % ридов, с одной ошибкой — 35 · 0,01 · 0,99³⁴ ≈ 25 %. Все 1908 легших ридов легли на то место, откуда их взяли. На поиск ушло в среднем 287 вызовов rank на рид, считая обе цепи и выдачу мест; простой проход по геному сравнивал бы рид с 97 004 местами.

Там же — история с HindIII. В эталоне слово AAGCTT встречается шесть раз, и куски получаются длиной 23 130, 2027, 2322, 9416, 564, 6682 и 4361 букв. Достаточно поменять одну букву, C на T в позиции 37 589, — и мест семь, а кусок 6682 распадается на 125 и 6557. Это ровно восемь кусков лабораторного маркера. С EcoRI и BamHI расхождения нет: шесть кусков EcoRI, 21 226, 4878, 5643, 7421, 5804 и 3530, — та же лесенка, что получают в лабораториях.

Код

Перестановку листинг строит через суффиксный массив удвоением, как в статье о нём: ранги пар упакованы в одно число, чтобы сортировать целые. Индекс — строка L, таблица C, счётчики через каждые 64 строки и позиции, кратные 32. Функция interval — обратный поиск, locate — шаги LF до отметки, one_mismatch — поиск с одной заменой.

Обратите внимание, что после построения программа больше не обращается к самому тексту: и поиск, и места, и проверка ридов идут по L, C, счётчикам и выборке. В Python счёт буквы внутри блока делает str.count — быстрый цикл внутри интерпретатора; в Go — обычный цикл по 64 байтам.

Сколько это стоит

Для ДНК σ = 4, и L занимает два бита на букву. В геноме человека около трёх миллиардов букв: суффиксный массив из 4-байтовых номеров занял бы около 12 ГБ, а программа Bowtie в 2009 году прикладывала к геному человека больше 25 миллионов ридов в час на одном ядре процессора и занимала около 1,3 ГБ памяти.

Построение — самая тяжёлая часть, но его делают один раз на геном. Удвоение из листинга для λ занимает доли секунды; для генома человека нужны линейные алгоритмы и построение по кускам, которое не держит весь суффиксный массив в памяти.

Варианты

Сжатие. Серии одинаковых букв в L хорошо сжимаются: bzip2 переставляет блоки до 900 КБ, перекодирует буквы «сколько разных было с прошлого раза» и сжимает результат кодом Хаффмана, — об этом глава курса «Конкурс упаковки». Сжатый индекс хранит L уже сжатой и умеет считать rank, не разжимая её целиком.

r-index. Если текст — тысячи почти одинаковых геномов одного вида, число серий r в L растёт медленно: новый геном почти не добавляет новых серий. В 2018 году Трэвис Гейги, Гонсало Наварро и Никола Прецца показали, как искать и находить места, храня O(r) чисел, а не O(n). У одного генома фага λ, как мы видели, серий почти столько же, сколько букв, — выигрыш появляется только на коллекциях.

Двунаправленный индекс. Обратный поиск продлевает образец только влево. Если построить ещё и преобразование перевёрнутого текста и согласовывать отрезки в обоих, образец можно продлевать в обе стороны. Так ищут с несколькими ошибками: рид режут на части, одна из которых наверняка без ошибок, ищут её точно и продлевают в стороны, тратя ошибки только там.

Выравнивание с пропусками. Риды бывают со вставками и пропусками букв, и одной заменой их не найти. Программы вроде BWA-MEM находят точные куски — затравки — обратным поиском, а потом выравнивают рид вокруг них динамикой, как в статье о выравнивании последовательностей.

История

Дэвид Уилер придумал перестановку в 1983 году, но не опубликовал. В 1994 году она вышла в отчёте Майкла Барроуза и Уилера «A block-sorting lossless data compression algorithm» исследовательского центра фирмы DEC как основа сжатия, и в 1996 году Джулиан Сьюард выпустил на ней bzip2. Указателем перестановку сделали Паоло Ферраджина и Джованни Манцини: в 2000 году они показали, что по сжатой перестановке можно искать за время, зависящее только от длины образца. Свою структуру они назвали FM-индексом — «Full-text index in Minute space», полнотекстовый указатель крохотного размера.

В биологию FM-индекс пришёл, когда машины стали выдавать миллионы коротких ридов за запуск. В 2009 году вышли Bowtie Бена Лэнгмида, Коула Трапнелла, Михая Попа и Стивена Зальцберга — с поиском с возвратом, который учитывает качество каждой буквы рида, — и BWA Хэна Ли и Ричарда Дурбина, Burrows–Wheeler Aligner. Их новые версии, Bowtie 2 и BWA-MEM, работают и сегодня, а r-index несёт ту же идею на коллекции из тысяч геномов.

Задачи

В обеих задачах дана только перестановка — текста нет. В первой нужно посчитать много слов быстрее, чем текст успеешь восстановить и прочитать. Во второй — найти места ридов с одной ошибкой, как в листинге.