Картирование коротких прочтений () — биоинформатический метод анализа результатов секвенирования второго поколения, состоящий в определении позиций в референсном геноме или транскриптоме, откуда с наибольшей вероятностью могло быть получено каждое конкретное короткое прочтение. Обычно является первой стадией в обработке данных в случае, если известен геном исследуемого организма.
Метод Платформы секвенирования нового поколения позволяют эффективно секвенировать миллионы коротких последовательностей длиной 50-500 пар нуклеотидов. Для этого молекула ДНК или кДНК разбивается на множество коротких отрезков, которые далее секвенируются параллельно. После получения отсеквенированных последовательностей этих коротких отрезков (прочтений) из них нужно восстановить полный геном или набор последовательностей кДНК. Для этого нужно определить для каждого прочтения наиболее вероятную позицию в референсном геноме.
В отличие от сборки de novo, которая собирает прочтения вместе для реконструкции до этого неизвестного генома, многие современные проекты имеют «референсный геном» — уже известный близкий геном другого организма. Или это может быть набор референсных последовательностей. В таком случае для придания прочтению какого-то смысла должна быть определена его позиция в референсных данных. Этот процесс и называется картированием () или выравниванием. Картирование может иметь несколько разный вид для разных задач, например: для геномного картирования должны отсутствовать большие гэпы, в то время как для 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 ищет способы соединить их с помощью интрона.