[ЗАСТАВКА] Добрый день!
Несмотря на смену декораций, мы продолжаем работать с снипами, то есть мы
ищем вариабельные позиции, используя короткие риды и референсный геном.
Давайте посмотрим, на чем мы остановились.
Мы остановились на том, что будем сейчас использовать команду samtools mpileup.
Ей для работы нужно два файла:
один файл — это ваш геном, второй файл — это bam-файл отсортированный,
обязательно отсортированный.
Если он будет неотсортированный, у нас ничего не получится, правда, программа
выдаст ошибку о том, что работать не может, потому что файл неотсортированный.
Давайте посмотрим, есть ли у нас.
Да, у нас есть оба файла.
Вот у нас бактериальный геном, и вот — bam-файл.
О, как видите, у нас есть замечательная опечатка в бактериальном файле.
Я вам покажу, как можно с помощью этого, если кто-то не знает,
взять и переименовать легко.
Мы просто берем и убираем отсюда букву d.
Отлично.
Получилось теперь у нас,
что есть два файла корректно названных, все замечательно.
Проверяем, что работает ли у нас samtools mpileup.
Работает.
Кроме samtools mpileup нам потребуется и еще
его вторая половинка — bcftools, которая раньше называлась, эта команда,
view, теперь, после нескольких обновлений, она называется call.
Собственно, как вы видите, при обоих запусках выдается набор подсказок,
и значит, это корректно установленная программа.
Прошу отдельно обратить внимание, когда будете устанавливать, взять последнюю
версию samtools, потому что за последний где-то год они обновлялись несколько раз,
и поначалу было много ошибок, поэтому реально возможны трудности,
если вы возьмете не ту или версию, которая находится в разработке.
То, что я делаю сейчас, сделано на последней версии samtools, и, в общем-то,
работает.
Нам нужны следующие команды.
На самом деле, она у меня даже есть, так что я его просто вам покажу, вот.
Идем. У нас есть samtools mpileup,
дальше мы ему говорим параметром −u, что нам нужно,
чтобы было uncompress — небинарный выход.
И формат файла геномного в fasta, это означает буква f.
Указываем бактериальный геном, по которому работаем, и говорим, в каком...
какой bam-файл использовать.
Если у вас несколько bam-файлов, то их тоже можно использовать одновременно.
На самом деле, если вы будете использовать несколько bam-файлов, то и output будет
тоже сделан для нескольких сразу образцов, то есть, на самом деле, если вы хотите
искать снипы в популяции, то вам нужно всю популяцию анализировать одновременно.
Это можно сделать, указав просто серию bam-файлов через пробел,
либо дав файл с параметром −b и имя файла,
в котором будет просто указан путь к каждому из bam-файлов.
После этого мы используем специальный символ «пайпа»,
собственно pipeline, для объединения двух программ.
И отправляем все, что у нас выдает samtools mpileup в bcftools call.
Параметры −c означает, что мы будем использовать старый алгоритм samtools,
потому что они в новой версии выпустили еще один для многоаллельных
локусов и редких вариантов, но нам это сейчас не нужно,
у нас будет просто consensus calling.
И параметр v означает, что мы ищем только вариабельные позиции.
Это мы указываем, в каком формате нам нужен output файл, и просто я показываю,
что можно вот его сделать, например, архивированным, указав z.
На самом деле, если вы не хотите, чтобы у вас файл сразу архивировался на выходе,
можно указать другой формат, использовав, например, обычное v, это будет...
вот видите, здесь есть список, каких у нас бывает, он может быть compressed,
uncompressed, compressed VCF и так далее.
Вот можно, в общем-то, сразу сделать uncompressed VCF,
кому как больше нравится.
И давайте попробуем запустить.
Отлично, все работает.
Как вы видите, первое, что сделал samtools,
он построил индекс для fasta-файла, потому что мы его ему не дали.
На самом деле индексы строятся очень быстро, поэтому можете на это не обращать
внимание, это вот времени много не занимает.
Дальше, в зависимости от мощности вашего компьютера, вам может потребоваться от,
я не знаю, 3-х, 5-ти, может быть 10-ти минут, на то,
чтобы программа закончила работать.
Программа закончила работать, так что мы можем продолжить.
Давайте посмотрим, что у нас получилось.
И как мы видим, что получилось 2 новых файла.
Один файл с индексом, который построился сразу,
и еще один файл, который называется variants call vcf gz.
Это заархивированный файл с вариативными позициями.
И многие программы могут его использовать уже в таком варианте,
мы можем сами попробовать что-нибудь из него простое посчитать.
На самом деле, есть вот такого вида командная строка, которая, на самом деле,
опять же использует bcftools и параметр filter,
для того чтобы отфильтровать все снипы с качеством больше 20.
Качество снипов, на самом деле, это тот же самый freight scale,
то есть качество больше 20 будет означать, что у нас...
то есть качество вообще в 20 означает, что у нас ошибка 1 на 100,
качество в 30 — 1 на 1000, качество в 40 — 1 на 10000,
и, в общем-то, я думаю, вы уловили тенденцию.
Больше 20, на самом деле, не то чтобы очень высокое качество,
но для небольшого датасета вполне нормально.
Опять же, когда вы фильтруете снипы или любые другие инделы, то, в общем-то,
понимаете, что вместе со снипами низкого качества неверными,
вы также отфильтровываете зачастую верные снипы с низким качеством.
Так что здесь всегда нужно приходить к некоторому компромиссу, вот.
Кроме bcftools filter, есть еще одна команда bcftools stats,
которая собирает некоторые базовые статистики, и мы посмотрим соотношение
транзиции к трансверсии, сделаем это с помощью вот этой команды.
Сколько это у нас времени займет?
Практически нисколько.
Ну собственно, вот набор базовых параметров, если кому-то это интересно,
можно опять же посмотреть,
какие дополнительные параметры выдает bcftools stats,
программа достаточно удобная, потому что вам не придется все это считать вручную.
Но мало ли, вам вообще интересно, как выглядит vcf-файл?
Давайте попробуем.
Смотрите, мы его можем взять и просто открыть как обычный текстовый файл.
То есть он, когда не упакован, его вполне можно смотреть глазами.
Из чего же он состоит?
В нем есть некоторый заголовок, который содержит метадату,
то есть здесь множество информации о том, которая генерится в
процессе создания этого файла и описывает некоторые поля из файла.
Кроме того, у нас есть набор колонок: есть обязательные поля, такие как...
имя хромосомы здесь называется, но могут быть скэфолды, контиги,
все что угодно, по сути, это именно имя фрагмента, из которого у нас это делается.
Дальше идет координата в этом фрагменте, референсная,
то есть координата фрагмента, ID фрагмента.
Да, собственно, ID, как вы видите, у нас здесь стоит точками,
потому что у нас они не пронотированы.
Если вы занимаетесь, например, снипами человека или любыми другими вариативными
позициями, у которых есть специальные, там, rsid, то vcf-файл можно дополнительно
редактировать и специальными программами вставить просто эти сюда номера.
У нас этой информации нет, поэтому у нас ставят просто точки.
Дальше у нас идет референсный аллель, то есть это буква,
которая была в референсном геноме.
После этого идет альтернативный аллель, то есть это буква,
которая является консенсусом в ридах, за которым идет качество.
Дополнительный фильтр поля, оно опять же не всегда используется, то есть здесь
стоят точки, но некоторые инструменты могут дополнительно вставить некоторые
пояснения в поле «фильтр», если у вас стояло, например, качество выше 20, там,
вы хотели оставить или еще что-то, то в поле может появиться просто, там, pass or
failed, то есть отметить, что эта позиция прошла какие-то фильтры или не прошла.
Дальше идет поле info,
которое может сильно различаться от выбранных параметров.
То есть вот оно, поле info.
И в зависимости от того, с какими параметрами вы запускали samtools,
будут здесь различные вещи отображаться.
Но из интересного и, думаю, всем понятного,
вот можете обратить внимание на вот это значение.
Это покрытие, покрытие в этом месте.
То есть и мы видим, что оно бывает разное: у нас есть как места с покрытием
очень низким, так и с покрытием значительно более высоким.
И обычно фильтровать нужно не только по качеству, но и по покрытию.
Дальше, у нас есть последние два поля, они соотносятся с колонкой,
вот которая имеет имя нашего файла, то есть sorted_bacteria.bam.
И если бы у нас было несколько bam-файлов,
то количество колонок уходило еще в последних, было бы значительно больше.
Сколько бы файлов дали, столько бы и было,
потому что эта информация конкретно о нашем образце уже.
И из самого интересного, здесь у нас есть описание генотипа.
Вот.
Генотипы в vcf-файле отображаются тремя вариантами.
1/1 значит, что программа предполагает,
что здесь была альтернативная гомозигота, то есть обе буквы, то есть обе аллели,
которые он нашел, отличаются от референсной.
Но как мы видим, здесь покрытие вообще небольшое,
поэтому он просто нашел одну и считает, что других вариантов не было.
Вот. В других вариантах с покрытием,
где было оно выше, то, в общем-то, там тоже значит,
что все буквы в этом месте отличались от референсной, то есть в ридах, всех ридах,
выровнявшихся на это место, по такой координате было отличие.
Вот. Кроме того, есть обозначение 0/0,
которое говорит о том, что в этом месте у нас нету отличия от референса.
В нашем варианте мы таких вообще не можем встретить,
потому что у нас всего один образец, и он, если не отличается от референса,
то программа его просто не будет показывать.
Если бы у нас было несколько образцов, то вполне возможна ситуация,
когда у одного бы было, например, полное отличие от референса (ну вот 1/1),
а у другого было бы, например, 0/0, потому что он не отличается.
Кроме того, есть еще гетерозиготы.
Конечно, когда мы говорим о бактериях, гетерозигот вообще быть не должно,
потому что это крайне странно, но сейчас мы посмотрим, есть ли они у нас или нет.
Отображаются они обычно 0/1.
Вот.
Кроме того, бывают обозначения, когда не слэш используется,
а просто вертикальная черта, это значит, что информация фазирована.
Но такого, в общем-то, я думаю, если встретите, то не сразу,
то есть «фазирована» — это когда мы знаем, на какой именно аллели лежит,
то есть на каком именно стрэнде лежит какая буква.
Давайте, во-первых, посмотрим, сколько у нас всего гетерозигот,
а дальше я скажу, то есть еще раз, почему они могут образовываться.
Вот наши файлы, вот он, 6 мегабайт занимает наш vcf- файл, но,
на самом деле, он, в общем-то, небольшой.
Во-первых, сколько всего в нем строчек?
То есть...
это практически нам скажет о том, сколько всего вариативных позиций там минус
заголовок и, там, метаинформация, которая была в начале, но ее там немного.
Ну то есть около 40 тысяч.
Отлично.
Теперь с помощью команды grep, а лучше использовать даже fgrep,
в текущем варианте он будет быстрее, нам все равно крайне...
мы смотрим, сколько у нас есть гетерозигот.
Ой, не то я сделал.
И добавляя параметр −с, мы смотрим, сколько у нас гетерозигот.
Их 200, то есть у нас, по крайней мере
программа нашла 200 гетерозиготных позиций в бактериальном геноме.
Собственно, тут, я думаю, каждый может попытаться бы хотя бы придумать или
подумать, почему это могло случиться, но я самые очевидные варианты сейчас озвучу,
смотрите.
Рид мог просто выровняться неверно, то есть если у нас некоторые повторяющиеся
области в геноме есть, то он мог, ошибочно выровнявшись с одной области на другую,
сказать, что это у нас здесь гетерозигота.
Дальше.
У нас может быть, например, следующий вариант,
когда просто есть смесь штаммов, то есть у нас, мы секвенировали не одну бактерию,
а несколько, и поэтому мы тоже можем получить гетерозиготные позиции,
но на самом деле здесь их количество крайне маленькое, и я думаю,
что если посмотреть подробней (то есть эти все риды же можно найти, посмотреть о том,
как они выровнялись), скорее всего это все-таки misalignment опять же,
но можно просто их по их содержанию посмотреть, вот.
Таким образом мы с вами получили файл с набором вариантов.
Давайте еще раз посмотрим, какие мы использовали, то есть мы использовали...
ну вообще, начали мы использовать с fastaq-файлы, которые вы почистили,
использовали геном, на который с помощью iii выровняли fastaq-файл, то есть
короткие последовательности, сделали sam-файл, который можно смотреть глазами,
и это просто набор координат для каждого из ридов.
После этого мы его перевели в бинарную форму, отсортировали и передали команде...
программе samtools для работы, вот, собственно,
которая в итоге сделала нам vcf-файл.
Для работы с vcf-файлом тоже есть инструмент, который, как ни удивительно,
называется vcftools, с помощью которого вы можете дополнительно фильтровать,
например, по качеству (если вас не устраивают,
там, фильтры, которые есть в bcftools), по покрытию или по регионам, в общем,
как вам удобно, используйте что удобней.
Но опять же обращаю внимание, samtools — не единственный вариант из тех,
что есть для текущей задачи.
Но он по крайней мере один из самых популярных.
Те ответы, которые вы будете получать, используя более-менее хорошо известный
инструмент, они будут пересекаться, но совпадать, конечно,
будут далеко не всегда, то есть именно полностью совпадать.
Но опять же в позициях, которые хорошо покрыты, с хорошим качеством, я думаю,
никаких проблем возникать не будет.
Вот. Собственно, мы подошли к тому,
чему хотели — мы получили файл с набором вариативных позиций.
И в дальнейшем мы попытаемся посмотреть, что еще можно с этим файлом сделать,
потому что это они так просто разбросаны по геному, а, например, мы хотели бы
знать, в каких именно генах они лежат или какие гены у нас самые вариабельные.
Но для того чтобы начать это делать, нам нужно на самом деле гены в геноме найти.
[ЗАСТАВКА]
