Курс «Введение в биоинформатику» адресован тем, кто хочет получить расширенное представление о том, что такое биоинформатика и как она помогает биологам и медикам в их работе. The course is aimed at those who would like to have a better idea of what bioinformatics is and how it helps biologists and medical scientists in research and clinical work.
Форматы SAM, BAM. Сортировка. (SAM, BAM Alignment Format. Sorting.)
Loading...
來自 Saint Petersburg State University 的課程
Введение в биоинформатику (Introduction to Bioinformatics)
45 個評分
從本節課中
Неделя 5 - Выравнивание коротких фрагментов (Short read alignment)
Пятая неделя нашего курса посвящена выравниванию коротких фрагментов на референс - подходу, который используется при работе с уже изученными геномами. В лекциях разбираются наиболее популярные программы для решения таких задач, а также основные алгоритмы для работы со строками. <br>Fifth week of our course is devoted to short read alignment to the reference genome, which is used when we work with the already studied organisms. In lectures we will talk about the most popular software for solving such problems, as well as basic algorithms for working with strings.
與講師見面
Alla L LapidusProfessor, Department of Cytology and Histology
Centre For Algorithmic Biotechnology
Николай Вяххи
Павел ДобрынинМладший научный сотрудник
Центр геномной биоинформатики им. Ф.Г.Добржанского
Михаил РайкоПостдок
Центр геномной биоинформатики им. Ф.Г.Добржанского
Екатерина ЧерняеваПостдок
Центр геномной биоинформатики им. Ф.Г.Добржанского
[ЗАСТАВКА] Для того чтобы у
нас все работало, нам нужно быть предельно внимательными и аккуратными.
Давайте начнем писать.
Первое.
Нам нужно имя программы, которую мы используем.
bowtie2.
Проверяем, что она работает.
Она работает.
После этого нам нужен параметр, о котором я упоминал,
это количество используемых ядер, −p.
И я хочу, чтобы сделалось все побыстрее, я поставлю 20.
После этого указываем индекс.
−x.
И индекс.
Мы называли его, начинали с bt.
Вот. Мы помним, что там получилось 6 файлов,
но на самом деле нам нужен не конкретный файл указать, а просто префикс, до точки.
То есть, вот bt2 bacteria.
Все. После этого, поскольку риды парные,
мы указываем ему левый и правый рид.
Первый будет −1 параметр просто и второй будет −2.
И также все замечательно.
−2.
И последний параметр −s, это у нас будет...
как мы его назовем: bacteria.sam.
Вот.
Запускаем.
В норме, если все работает хорошо и правильно, и вы запускаете на своей
локальной машине, он у вас должен повиснуть в ожидании и ничего не писать.
У меня он выдал предупреждение, но у вас, скорее всего, такого не будет, оно вообще
не связано даже с этой программой, так что на это не обращайте внимание.
Это должно закончиться предельно быстро, то есть у нас сейчас, я думаю, буквально,
может, даже минуты нету, и у нас получится sam-файл.
Также, кроме того, программа должна выдать некоторый набор статистик.
После выравнивания она скажет, сколько ридов у нас выровнялось,
сколько ридов выровнялось всего, сколько из них выровнялось совместно,
то есть которые действительно...
О! Все готово.
И давайте посмотрим.
Первое.
Вот. Пройдемся по тому, что получилось.
У нас, во-первых, все риды парные, но мы это уже проверяли, когда посмотрели, что у
нас количество ридов в общем-то совпадает в обоих файлах, то есть все парное.
Дальше.
Вот смотрите.
Первая, казалось бы, чудовищная проблема.
У нас 60 % из всех ридов не выровнялись совместно.
То есть они вроде как парные, но в то же время как пары выровняться не захотели.
Странно.
Но. Хорошая новость.
40 % ридов практически выровнялись, и только один раз.
То есть, как минимум 40 % всей нашей библиотеки это именно то, что нам нужно.
Есть некоторое количество повторов среди тех, которые
выровнялись вместе в несколько раз, то есть, ну, крайне небольшое количество.
Видите, он просто говорит, что больше одного раза у нас выровнялось 0,5 %.
Ничего страшного.
Хорошая новость, что из тех, которые не выровнялись как парные,
они выровнялись по одиночке еще 60 % из них.
То есть, в общем-то, у нас не выровнялось не так много,
как могло показаться сначала, но какое-то количество,
то есть, у нас библиотека оказалась не то чтобы прямо очень парной.
То есть, скорее всего, в приготовлении могла произойти ошибка,
и у нас образовались некоторые химеры, либо у нас референсный геном,
на который мы выравнивали, тоже оказался не слишком подходящий, возможно,
здесь различия больше, чем ожидались.
Вот. В идеале мы могли бы ожидать в некоторых
случаях 98—99 % ридов могут выровняться идеально, хорошо и прекрасно.
Но бывает всякое.
Бывает и 40 %, может быть и 15 % в некоторых случаях.
Самое хорошее, что в итоге у нас все получилось.
А именно, у нас получился сам файл.
Давайте на него глянем.
Вот он у нас.
И сразу давайте обратим внимание на объемы.
У нас сам файл по размеру практически совпадает сумме двух файлов с ридами.
То есть, он может быть чуть больше, чуть меньше, зависит от того,
отсортировано или нет.
Мы это посмотрим.
Но просто будьте готовы к тому, что, например,
если у вас был геном бактериальный, пускай всего там несколько мегабайт, еще пару,
может быть, иногда риды могут занимать всего несколько гигабайт, иногда это могут
быть десятки или сотни гигабайт, то после того, как вы сделаете сам файл, вам нужно
еще столько же места, сколько было, просто для того чтобы хранить этот файл.
Внутри он в общем-то...
Есть отдельная на него спецификация, но мы посмотрим просто некоторые поля,
которые могут быть интересны.
Здесь, на самом деле, можно сказать, некоторые усложнения из того,
что было в fast equ.
Единственное могу сказать, что с помощью небольшой хитрости и минимальных
навыков программирования мы sam-файл можем превратить опять в fast equ.
Потому что у нас здесь есть и последовательность нуклеотидов в риде,
и качество, ассоциированное с этой последовательностью.
Вот.
Кроме того, у нас есть заголовок,
а именно, откуда у нас были взяты риды.
У нас есть дополнительный флаг на самом деле, который здесь различные
просто цифровые кодировки, можно, есть отдельный список, какие они бывают.
И здесь закодировано, что у нас рид выровнялся уникально,
не уникально, парно, непарно, в разные стороны, и так далее.
У нас эта информация в файле есть, но мы ей пока пользоваться не будем.
Дальше.
У нас есть заголовок того, у нас поскольку файл был весь геномный в одном фрагменте,
это fast-файл, то здесь у всех ридов будет один заголовок.
Если бы у вас было бы там multifast, то здесь могли бы быть разные заголовки.
И координата куда выровнялась?
Там есть некоторый набор дополнительных параметров, который говорит о том,
насколько хорошо выровнялось.
На самом деле вот эта следующая колонка цифр — это именно качество выравнивания.
То есть, на самом деле, вы сам файл потом можете тоже отфильтровать, например,
дополнительно по качеству, если вы считаете,
что у вас оно должно быть не ниже какого-то уровня.
Качество, на самом деле, очень похожее на то,
с чем вы уже сталкивались при оценке качества ридов.
Вот.
И вы получили sam-файл.
И для того чтобы работать с ним, нам понадобится уже другая программа,
которая называется samtools.
Для дальнейшей работы нам потребуется новый пакет, а именно, набор утилит,
samtools который называется,
разработанный специально для работы с файлами формата sam.
Давайте проверим, что у вас все работает хорошо.
Для этого можно просто набрать samtools.
Вот. Дает некоторую инструкцию, по идее, о том,
какие внутри него есть подпрограммы и что можно использовать.
И просто покажу хороший пример, что можно сделать с вновь полученным sam-файлом.
На самом деле sam-файл большой, и его так хранить не очень удобно.
И лучше его перевести в другой родственный формат, который называется bam,
от «b», потому что binary, и его уже нельзя будет смотреть глазами,
потому что это будет бинарный файл, но занимать места он будет меньше.
Давайте посмотрим.
Для этого нам потребуется команда samtools view и 2 аргумента,
которые нужны, это −S большое, говорит о том,
что файл на входе будет как sam, а b — чтобы его на выходе сделали bam-файлом.
Указываем имя файла
и какой файл мы хотим записать.
По дефолту программа будет выводить его в стандартный output,
то есть она будет на экран выводить, а это нам не нужно,
поэтому мы ставим значок записи в файл и записываем в bam.
Этим можно ограничиться, но есть еще один маленький момент,
который можно использовать.
Дело в том, что samtools view команда, точнее,
набор команд тут будет внутри, которую можно использовать, не
ограничивается только перегоном из sam в bam, там есть еще несколько разных вещей,
которые вы можете посмотреть, набрав просто help.
Но одна из очень полезных это, например, F большое и 4.
F — это «отфильтровать»,
а 4 — это цифровой код сигнала в sam-файле,
означающий, что рид никуда не выровнялся.
В общем-то нам для дальнейшей работы совершенно не важны те риды,
которые не выровнялись на геномы, а они у нас в sam-файле есть и занимают место.
Таким образом, мы при конвертации из одного формата в другой можем вполне
успешно решить еще одну задачу и выкинуть все ненужное.
Но если вы хотите наоборот посмотреть, что же это у вас туда не выровнялось,
потому что причины, по которым у вас не выровнялись риды на геном,
который вы исследуете, могут быть совершенно разные.
Это может быть другой штам, и вы хотите узнать, что в нем другого, это может быть
примесь другой бактерии, и вы опять же хотите узнать, что же там было такое.
То есть, в общем-то, это один из вариантов того, как можно отсеять риды,
которые вам не нужны или, наоборот, нужны.
Поэтому ставим −F большое, если поставить f малое, то, наоборот,
останутся только те, которые не выровнялись.
Если F большое, то их, наоборот, уберут.
Вот. Цифровых кодов там совершенно много,
и опять же их можно найти просто в Интернете, в инструкции к sam-формату.
Теперь запускаем.
Программа начала работать и проработает какое-то время, но это должно быть
достаточно быстро, опять же все зависит от мощности отдельно взятого компьютера.
Из минусов могу сказать,
что в этом случае как раз samtools не очень хорошо распараллеливается.
И придется ждать, пока он закончит.
Времени это, правда, особо много занять не должно,
потому что у нас бактериальный геном, и файлы не очень большие.
Вот.
Кроме того, в samtools я буду использовать только несколько в текущей
демонстрации вещей, которые можно сделать (вот, кстати, он закончил работать).
Но на самом деле samtools можно с помощью этой программы считать глубину покрытия,
степень покрытия и объединять несколько sam-файлов.
Наоборот, если у вас есть, например, был большой файл с кучей хромосом,
а вы хотите достать только одну хромосому, это тоже можно сделать с помощью samtools.
Вот. Ну, у нас все закончило работать.
Давайте посмотрим,что у нас получилось.
[ПЕЧАТАЕТ] Мы видим,
что у нас есть наш sam-файл, с которым мы начали работать.
И вот получился bam-файл.
Как видите, он значительно меньше.
То есть, в среднем можно ожидать, что он уменьшаться будет в 4—5 раз.
На самом деле ничего удивительного, потому что внутри,
по сути, формата bam это тот же самый sam, который
просто был дополнительно переупакован обычным архиватором типа kzip,
который вы могли бы ожидать на любой линуксовой машине.
Вот. Но с этим bam-файлом пока еще рано
что-либо делать еще, потому что он не сортированный.
Риды в нем расположены в том порядке, в котором они были у вас в файле с ридами.
То есть каждая пара идет, вот начиная с той,
которая была верхней в файле с ридами, та и оказалась сверху.
Это очень неудобно, и в дальнейшей работе нам такое не годится.
Для того чтобы отсортировать, есть специальная команда,
которая называется samtools sort.
И опять же инструкция к ней крайне простая.
Нам всего лишь нужно сказать, что мы сортируем.
Это у нас будет bacteria.bam и префикс того файла,
который получится после сортировки.
Это будет sorted bacteria.
Если вы впишете sorted bacteria.bam, то он потом пришьет еще один .bam.
Поэтому в общем, ничего смертельного, можно потом переименовать,
но здесь просится указывать именно префикс.
Вот.
И давайте запустим.
Тут будет отдельный момент, связанный с сортировкой.
Во-первых, это позволит дальше работать и является вообще необходимым
шагом для всего остального.
После сортировки у нас все риды будут расположены относительно координат в
геноме, то есть сначала будут идти риды, которые выровнялись с координатой,
грубо говоря, 1.
Потом следующие риды, которые выровнялись на координату 2, 3 и так далее.
Но уже не будет соответствовать тому порядку, который был в исходном файле,
но зато очень удобно для работы в дальнейшем.
Вот. Уже закончил работать.
Опять же на больших геномах это может занимать,
или если у вас намного больше объем данных, значительно больше времени.
Прямо намного больше времени.
Но тут это все пока достаточно быстро.
В процессе работы опять же программа будет создавать несколько временных файлов,
так что не пугайтесь, если там будет, например,
называться что-нибудь sorted bacteria 001, 002 или 003 и так далее,
она их потом объединяет все в один файл и ненужные файлы удаляет.
Вот.
Давайте еще раз посмотрим, что у нас получилось.
Какие bam файлы у нас теперь есть.
Момент.
Вот так.
И удивительное дело!
Сортированный файл стал еще меньше, причем значительно меньше, чем несортированный.
Тут опять же пугаться не стоит, это нормально.
После того, как файл стал сортированный,
возможность его эффективно заархивировать только увеличилась.
Именно поэтому он стал меньше.
На самом деле потому что некоторые начинают бояться, что что-то стало,
они куда-то потеряли данные и так далее.
Нет. Это абсолютно нормально, и сортированный
файл всегда будет практически занимать меньше места, чем не сортированный.
Так что вы все сделали правильно, если у вас все именно так и выглядит.
Опять же, для примера, как можно проверить?
То есть, если вы попытаетесь открыть этот bam-файл,
то вы ничего вразумительного не увидите,
потому что этот формат не предназначен для чтения людьми.
Он только для компьютера.
Вот. Но вы можете все равно, во-первых,
его обрабатывать с помощью других программ, а к тому же, если вы хотите
узнать, что же у вас там было, то вы можете с помощью той же команды samtools
view посмотреть, например, заголовки sequence, которые там есть.
То есть, это именно те заголовки,
которые были у multifast и на которые выравнивались.
То есть, это будут, если это был большой геном, например,
эукариотически это могут быть имена хромосом или скефолдов,там,
или имена бактериальных хромосом, если их было несколько.
Вот.
А! И как видите, это заголовок bam-файла.
Здесь описаны программы, с помощью которых он был получен исходно, то есть сам файл,
версии и имя sequence,
на который выравнивалось, плюс его размер.
Вот. На текущий момент у нас все получается
корректно, надеюсь, у вас тоже.
Следующим шагом будет как раз уже, поскольку мы в
текущем задании будем пытаться получать снипы или SNP,
или однонуклеотидные полиморфизмы, или ОНП, как их иногда называют по-русски, то
дальше как раз будет следующая команда, с помощью которой это и можно будет сделать.
[ЗАСТАВКА]
