Картирование коротких прочтений ( англ. Short-Read Sequence Alignment, Short-Read Sequence Mapping ) — биоинформатический метод анализа результатов , состоящий в определении позиций в референсном геноме или транскриптоме, откуда с наибольшей вероятностью могло быть получено каждое конкретное короткое прочтение. Обычно является первой стадией в обработке данных в случае, если известен геном исследуемого организма.

Метод

Платформы секвенирования нового поколения позволяют эффективно секвенировать миллионы коротких последовательностей длиной 50-500 пар нуклеотидов. Для этого молекула ДНК или кДНК разбивается на множество коротких отрезков, которые далее секвенируются параллельно. После получения отсеквенированных последовательностей этих коротких отрезков (прочтений) из них нужно восстановить полный геном или набор последовательностей кДНК. Для этого нужно определить для каждого прочтения наиболее вероятную позицию в референсном геноме.

В отличие от сборки de novo , которая собирает прочтения вместе для реконструкции до этого неизвестного генома, многие современные проекты имеют «референсный геном» — уже известный близкий геном другого организма. Или это может быть набор референсных последовательностей. В таком случае для придания прочтению какого-то смысла должна быть определена его позиция в референсных данных. Этот процесс и называется картированием ( англ. mapping ) или выравниванием. Картирование может иметь несколько разный вид для разных задач, например: для геномного картирования должны отсутствовать большие гэпы, в то время как для RNA-seq они допустимы из-за наличия сплайсинга. В общем и целом, задачи картирования не изменились с прошлого поколения секвенаторов, но программы, разработанные для предыдущего поколения, не предназначены для работы с возросшими объёмами данных, а также плохо работают с прочтениями короткой длины.

Проблемы

Вариабельность генома и ошибки секвенирования

Главная проблема состоит в том, что исследуемый геном отличается от референсного из-за вариабельности генома (например, из-за SNP или инделов), а также из-за ошибок секвенирования. Из-за этого выравнивание прочтения и его «правильного» положения в референсном геноме может иметь больше различий между собой, чем выравнивание этого участка на какое-либо другое место в референсном геноме, поэтому программы картирования должны уметь находить неточные соответствия. Для этого используются различные эвристики. При наложении на это возможности сплайсинга в случае с RNA-seq проблема становится ещё сложнее.

Повторяющиеся элементы

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

Вычислительная проблема

Если референсный геном велик и сделаны миллиарды прочтений, то серьёзной проблемой является время картирования. Выравнивание всегда было крайне ресурсоемкой задачей, но в данном случае проблема выходит на качественно новый уровень, требуя от алгоритмов экстраординарной рациональности и эффективности использования процессорного времени и памяти.

Подходы

Существуют два основных подхода к решению этих задач: с использованием хеш-таблиц и с использованием суффиксных деревьев или массивов.

Основы подхода с использованием хеширования

Процесс поиска с помощью хеширования во много раз быстрее и менее затратный, чем классическое выравнивание с помощью динамического программирования для алгоритма Смита — Ватермана .

В данном подходе используется хеш -функция, которая позволяет трансформировать строку в ключ, который может быть использован для быстрого поиска. Наиболее простым способом было бы разбиение генома на слова, совпадающие по длине с прочтением, но этот подход не работает, так как длинные слова будут уникальны с большей вероятностью и их хранение занимало бы слишком много места в памяти. Вместо этого используются хеширование более коротких участков, которые встречаются гораздо чаще. После того, как с помощью хеш-функции будут получены подходящие позиции, можно попробовать картировать оставшуюся часть прочтения на геном. Также подход разделения прочтения на несколько частей позволяет заложить в алгоритм возможность замен. Так в программе Maq прочтение делится на 4 части, называемые сидами (короткие участки с точными совпадениями). Если прочтение идеально ложится на референс, то идеально ложатся и все 4 сида. Если в прочтении есть одна замена (mismatch), которая, вероятно появилась из-за наличия SNP или ошибок секвенирования, то она попадает в один из сидов, в то время как остальные 3 будут по-прежнему идеально картироваться. Аналогично, при двух заменах минимум два сида будут идеально картироваться. Таким образом, картировав с помощью индексирования все возможные пары прочтений (их будет 6, так как сиды должны идти в определённом порядке), мы получим ограниченный набор мест в геноме, куда может быть картировано все прочтение, после чего возможно использовать уже обычное выравнивание для определения того, какая из позиций подходит лучше всего (смотри изображение 1а). Схожим образом работают SOAP, RMAP и SeqMap.

Модификацией данного подхода может служить идея рассматривать все k-меры прочтения, вместо неперекрывающихся участков определённой длины. Так для прочтения ACGT их будет 3: AC, CG, GT. Подобным образом работают SHRiMP2 и RazerS. Это увеличит чувствительность, но увеличит и затраты времени.

Вся эта информация занимает достаточно много места в памяти. Для уменьшения объёма потребляемой памяти большинство программ используют двухбитное кодирование нуклеотидов (A 00, C 01, G 10, T 11), но это не позволяет использовать неоднозначный нуклеотид N, который может присутствовать как в прочтениях, так и в референсном геноме. Программы используют разные подходы, чтобы обойти это. Так BWA заменяет N на случайный нуклеотид, а SOAP2 заменяет все N на G.

Для ускорения работы алгоритмов и обхождения ошибок можно использовать различные улучшения. Например, использовать индифферентные позиции (обозначим их X). Так сид ACGxACG будет выравнен и на ACGAACG, и на ACGCACG, что увеличивает чувствительность картирования (правда, увеличивает время работы, так как на каждое прочтение будет больше находок). Это используется в таких программах как Zoom, BFAST, GASSST, SHRiMP2, PerM.

Большую часть времени алгоритмы тратят не на поиск сидов, а на проверку их окружения. Большинство программ используют Алгоритм Нидлмана — Вунша или его модификации. Другие, например GASSST, добавляют промежуточный шаг, заключающийся в измерении Эйлерова расстояния, которое учитывает просто количество одинаковых букв. К примеру, если мы пытаемся картировать прочтение, содержащее 5 G на участок, содержащий только 1 G, то мы будем иметь как минимум 4 замены. Этот подход позволяет быстро отсеять неподходящие участки и применять более аккуратные (но и затратные) подходы только к перспективным регионам.

Возможно хешировать не геном и искать в нём прочтения, а хешировать прочтения и искать в них участки генома такой же длины. В ранних версиях Maq, Rmap и SeqMap использовался такой подход, но с тех пор количество прочтений в одном эксперименте сильно возросло и этот подход перестал быть рациональным.

Основы подхода с использованием суффиксов

Алгоритмы, основанные на хешировании, плохо справляются с повторами, так как сильно возрастает количество сидов, которые приходится проверять. Для решения этого используются алгоритмы, основанные на суффиксных деревьях и суффиксных массивах . Преимущество данного подхода, в частности, заключается в том, что повторы не увеличивают время работы алгоритма, так как повторяющиеся участки «схлопываются» в суффиксном дереве. В чистом виде этот подход будет работать крайне быстро при условии отсутствия ошибок и замен (например, его использует программа MPscan).

Суффиксный массив используется для большей экономии памяти. В общем и целом, поиск по суффиксному массиву принципиально не отличается от работы с суффиксным деревом, но требует несколько более сложного подхода. Часто используется трансформация Барроуза — Уилера . После всех преобразований размер такого суффиксного массива сопоставим с размером исходного генома. Так для всего человеческого генома суффиксный массив, создаваемый программой bowtie, занимает 2 гигабайта. Для сравнения, база данных индексов, основанных на хешах (как та, что используется в программе Maq), занимает около 50 гигабайт памяти.

Для поиска слов используется .

В упрощенном виде процесс заключается в следующем. Прочтения выравниваются по одному нуклеотиду к геному, преобразованному по Барроузу — Уилеру. Каждый выравненный нуклеотид позволяет программе сузить количество мест, куда может ложиться все прочтение. Если программа не может подобрать позицию, куда прочтение ложится идеально, она возвращается к предыдущему шагу, делает замену и продолжает поиск. Таким образом работает, например, SHRiMP. С другой стороны, если количество разрешенных ошибок велико, то это начинает замедлять алгоритм. Чуть более интересная модификация используется в программе BWA — она сначала выравнивает левую и правые части прочтения, предполагая, что хотя бы в одном из этих двух регионов будет меньшее количество ошибок (что основывается на факте, что 5'-конец обычно секвенирован лучше).

Сравнение подходов

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

RNA-seq

В данном случае в рассмотрение следует включать возможность сплайсинга. В основном используется знание об уже известных экзонах и интронах, что позволяет создать базу данных, состоящую из экзон-экзонных объединений, и уже на ней реализовывать обычные алгоритмы, такие как bowtie или maq. Очевидным образом, этот подход не учитывает неизвестные ранее интроны, поэтому он может быть совмещен с машинным обучением для предсказания неизвестных сплайсингов.

Другой подход может не использовать аннотацию вообще. В режиме работы без аннотации TopHat на начальном этапе определяет потенциальные экзоны, строит по этой информации базу потенциальных сайтов сплайсинга и затем производит картирование прочтений с её использованием. Потенциальные экзоны определяются по местам смежных RNA-seq прочтений при выравнивании их на геном. Поскольку многие экзоны короче, чем генерируемые современными секвенаторами прочтения, они не будут обнаружены во время начальной фазы картирования. TopHat решает эту проблему в основном с помощью разбиения прочтений на более короткие фрагменты и их независимого картирования.

TopHat использует два основных метода при определении потенциальных сайтов сплайсинга. Первый и основной фактор в пользу потенциального сайта сплайсинга — когда два сегмента с одного прочтения (для прочтений длины 45 пар нуклеотидов и более) картируются на определённом расстоянии друг от друга или когда не удаётся картировать внутренний сегмент прочтения. Другой фактор — появление пар «островов покрытия», которые являются местами, где много RNA-seq прочтений картируются рядом. Соседние острова часто вырезаются вместе и находятся рядом в транскрипте, поэтому TopHat ищет способы соединить их с помощью интрона.

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

Программы

BFAST

Аббревиатура от англ. Blat-like Fast Accurate Search Tool . Разработчики программы сделали акцент на чувствительность программы к ошибкам, SNP и инделам, позволяя выбирать баланс между скоростью работы и точностью.

Поддерживает секвенирование спаренных концов . Использует для финального выравнивания алгоритм с поддержкой гэпов и определением небольших инделов . Умеет работать в параллельном режиме на кластере. Есть версия программы bfast+bwa. Поддерживает форматы данных Illumina, ABI SOLiD, 454, Helicos.

BLAT

BLAST-подобное средство выравнивания. Оптимизирован для скорости за счёт использования индекса непересекающихся K-фрагментов, помещающегося в оперативную память компьютера .

Допускает одну замену в каждом совпадении.

Bowtie

Использует алгоритм Барроуза—Уилера для индексирования. Программа оптимизирована по скорости и расходу памяти, умеет использовать несколько ядер процессора. Заявлена скорость превышающая в 35 раз скорость MAQ и в 300 раз — скорость SOAP в одинаковых условиях.

Допускает несовпадения.

На базе Bowtie построена программа TopHat для выравнивания прочтений RNA-seq.

BWA

BWA (выравнивание биологических последовательностей) — это комплект из трёх программ: BWA-backtrack, BWA-SW and BWA-MEM. BWA-backtrack работает с прочтениями от Illumina длиной до 100 пар нуклеотидов, BWA-SW и BWA-MEM могут работать с более длинными последовательностями от 70 до 1 миллиона пар нуклеотидов. BWA-MEM является последней версией, более качественной и более точной.

BWA-SW и BWA-MEM умеют находить химерные прочтения.

BWA-SW использует преобразование Барроуза—Уилера вместе в выравниванием Смита—Ватермана. Умеет работать с длинными последовательностями, при этом точнее и быстрее, чем BLAT.

ELAND

Расшифровывается как Efficient Local Alignment of Nucleotide Data. Разработка компании Solexa, приобретённой затем Illumina.

Использует paired-end-прочтения, умеет определять структурные варианты. Умеет работать только с последовательностями длиной до 32 пар нуклеотидов, допускает до двух различий, но не может работать с инделами.

MAQ

Производит выравнивание без гэпов. Для single-end-прочтений может найти до 2 или 3 несовпадений в зависимости от параметров запуска. Для paired-end-прочтений допускает до 1 несовпадения.

Строит консенсус на основе статистической модели.

Novoalign

Выравнивает single-end- и paired-end-прочтения от Illumina, paired-end от 454. Может определять выравнивания с гэпами или несовпадениями. Может сообщать о множественном выравнивании.

SHRiMP

SHRiMP2 делает акцент на точности, позволяя выравнивать прочтения с полиморфизмами или ошибками секвенирования.

Использует алгоритм Смита-Ватермана. Версия 1 индексировала прочтения, версия 2 индексирует геном, за счёт чего достигает большей скорости.

Поддерживает прочтения Illumina/Solexa, Roche/454 и AB/SOLiD, поддерживает параллельные вычисления.

SOAP

Умеет производить выравнивание single-read и pair-end фрагментов. Умеет работать с интронами.

Алгоритм использует индекс 2way-BWT (2BWT) . Версия SOAP3 оптимизирована для работы с GPU и использует специальный индекс GPU-2BWT .

TopHat

Программа для выравнивания прочтений RNA-seq, построена на базе Bowtie.

Была создана для работы с прочтениями, производимыми Illumina Genome Analyzer, но успешно используется и с прочтениями, генерируемыми другими технологиями. Оптимизирована для прочтений длины 75 пар оснований и больше. Не позволяет смешивать вместе paired- и single-end-прочтения.

Сравнительные характеристики

Программа	Алгоритм	paired-end/single-end	Гэпы (интроны)	Инделы	Замены	Длина прочтений (п.н.)
BFAST	Смита—Ватермана. Есть версия с BWT	+/+	+		+
BLAT	K-меры (как BLAST)		+	1-2	1-2
Bowtie	Барроуза—Уилера	-/+			+	<=100 , 70-1M
BWA	Барроуза—Уилера + Смита—Ватермана	+/+	+
ELAND	Хеширование сидов?	+/?	-		<=2	до 32
MAQ	Хеширование сидов	+/+	-	+	2-3 для single-end, 1 для paired-end	<=63
Novoalign		+/+	+		+
SHRiMP	K-меры + Смита—Ватермана	+/+	+	+	+
SOAP	Барроуза—Уилера	+/+	+		<=2	<=60
TopHat	Барроуза—Уилера	+/+	+	+	<=2	>=75

См. также

• Высокоэффективное секвенирование
• Методы секвенирования нового поколения
• Полони-секвенирование

Примечания