И добро пожаловать на 6-ю неделю, на которой, мы сейчас с вами попробуем сделать аннотацию генома, используя две различных программы, которые ищут гены генома. Как вы знаете, у нас аннотация может быть генома, и может быть основана на гомологии. Генома программы, они работают зачастую значительно быстрее и потребуют от ваших компьютеров меньше ресурсов и вы сможете большим количеством времени поэкспериментировать и с большим количеством попыток поэкспериментировать сами, например, с различными параметрами. Я покажу вам, как можно воспользоваться программами «глиммер», которые используют претренированную модель для поиска и «джинмарк с», которой вообще, в общем-то она с самообучением, ей не нужно давать никаких моделей. Работать мы будем с тем же самым геномом, с которым мы работали на прошлой неделе и для которого у нас уже есть снипы. И, смотрим еще раз, давайте. Что у нас осталось с прошлой недели? С прошлой недели у нас осталось пять файлов, из которых нам, на самом деле, сейчас понадобится только файл с геномом и файл со снипами в конце. Давайте, создадим две новых директории. Одну мы назовем «глиммер». «Глиммер», а вторую мы назовем «джинмарк с». И начнем с «глиммер». Вам его надо будет установить самим. Я надеюсь, что к этому времени, у вас уже значительно больше опыта, хотя, понимаю, что зачастую это может быть очень трудно, разбираться с установкой всех этих программ. Но будьте готовы, что зачастую, эти самые большие трудности в работе могут быть связаны с тем, что просто не удается поставить или понять, по какой причине оно вообще не работает. «Глиммер» в этом смысле, программа достаточно простая. Одна трудность, кроме самого «глиммера», нужно поставить будет другую программу, которая называется ELF – это часть «глиммера», отвечающая за анализ мотивов. Вот, но в чем прелесть «глиммера» в том, что у него внутри уже есть некоторый набор скриптов, которые позволяют автоматизировать процесс аннотации. В том числе, есть скрипт, который называется «джитрирейтед» который практически сделает все сам. Давайте, проверим, что он работает, да. Он работает, правда, ругается, потому что, я ему ничего больше не дал. А от нас он хочет, смотрите что. Он от нас хочет геном и тег. То есть, то, как будут называться файлы, которые он произведет. Смотрим, папка у нас пока пустая. И давайте, запускаем. То есть, у нас в директории выше лежит бактериальный геном и мы это назовем, и мы это назовем, давайте, «глиммерджинс». И поехали. Как вы видите, он начал работать. Будет это он это делать не так долго, так что, давайте просто подождем, когда он закончит. Да, он может начать выдавать ошибки, с которыми придется разбираться отдельно, но обычно помогает чтение различных форумов. Даже если вы еще не сталкивались, отлично помогает «биостарс» «секансерс» и например, «старковерфол». Если у вас какие-то проблемы с компиляцией или еще чем-то. Количество проблем, на самом деле, может быть достаточно большим и по разной причине, так что, все их предугадать и ответить, может быть достаточно сложно, но если вы используете «Линукс», то наверное, у вас большей части проблем будет меньше. Если вы будете пытаться это устанавливать и компилировать на «МАК», то могут быть некоторые трудности, но как видите, ничего непреодолимого быть не должно. Так что, просто ждем окончания. Как видите, это не заняло много времени. Давайте посмотрим, что у нас получилось. Как всегда, мы используем ту же самую команду «олас». О, огромное количество файлов самых разных, но нам интересней всего, на самом деле, думаю, что будут файлы с координатами генов и вы можете догадаться, что они называются «глиммер», «джинс», «кордс». Вот, и смотрите, он в этот файл складывает открытые рамки считывания, которые он нашел и которые он предполагает, что гены. Также, он им считает... различный скоринг вводит для оценки того, что он сделал. В общем-то получилось все неплохо. Давайте посмотрим, сколько всего он смог найти таких файлов, ой, не файлов конечно, генов. Всего, он нашел 3, 200 с небольшим генов. Отлично. Собственно, на самом деле, смотрите, у нас уже фактически, то есть, один вариант аннотации есть. Хорошо, переходим к другой программе «джинмарк с». Устанавливать его даже проще, чем «глиммер», там просто меньше всего делать надо. И он уже полностью готов к работе. Смотрите, он опять же отозвался на мою команду и несмотря на обилие параметров, ему нужно не так много для работы. Ему нужно только «сиквенс» файл и все, в «фаста» формате. Прекрасно. На самом деле, именно это мы ему можем дать очень легко. И как обычно, запускаем. Программа закончила работать. Давайте посмотрим, что нам выдало «джинмарк». Файлов еще больше. Вот, но смотреть в этот раз стоит на бактерия «дженом фаста ист». Давайте заглянем внутрь. О, и еще один замечательный файл, в котором у нас есть координаты того, что нашлось. Собственно, у нас есть просто порядковый номер, у нас указан стренд, есть левый и правый, ну, собственно край, то есть начало... начало и конец, вот. Плюс к этому опять же у нас есть имя хромосомы, с которой это все было получено. Вот... э-э-э... давайте посмотрим общее количество того, что получилось. Поскольку они пронумерованы, нам даже не нужно пересчитывать, мы можем просто, мы можем просто его взять и показать хвост файла. О, удивительно. Получилось, на самом деле, приблизительно такое же количество с помощью двух разных программ, вот. Теперь, было бы очень интересно узнать, а насколько это все соотносится вообще с реальностью. Для этого, во-первых, на самом деле, вы могли бы взять координаты файла, аннотация оригинального, то есть, имена хромосомы есть? Есть Можно вылезти в «ин си би ай» и найти все эти гены. Или в любом другом месте найти аннотацию. Но, на самом деле, если посмотреть сколько всего, аннотированного генов в, в этом фрагменте, я имею в виду даже в официальной аннотации, то там как раз будет число, близкое к 3 100, там может 3 126. Вот, что уже очень интересно, то есть, у нас оба инструмента выдали достаточно точный результат. По крайней мере, он похож на то, что было найдено до этого. И теперь, на самом деле, очень интересно посмотреть, а насколько они вообще сходятся между собой, то есть, у нас нашлись какие-то фрагменты и интервалы а как проверить, что они нашли что-то похожее. Может это просто они подобные куски приблизительно в том же самом количестве. На самом деле, для этого вам нужно использовать будет следующий инструмент, который называется «интерсектбед», на самом деле, можете и не использовать этот инструмент. То есть, здесь крайне простая задача. Если кто-то умеет программировать, может вполне с этим справиться сам. То есть, нам нужно посчитать пересечение двух интервалов, то есть, набора интервалов и выдать какую-то статистику о том, насколько они пересеклись. Но, вот «интерсектбед» умеет это все делать сам. Как видите, он установлен и работает. И работает крайне просто, собственно, вот форматы файлов, с которыми он может работать - - «бед», «джфф» и «всф». «Всф» файл мы уже с вами встречали. «Джфф» файл – это стандартный файл для работы с генами. То есть, там есть его координаты, начало, то есть, обычно они там идут, ну еще, то есть, расписаны зоны. «Бед» файл тоже бывает разный, но «бед» файлы наиболее простые из всех этих вариантов. То есть, «джфф» файл, он сложнее по структуре, хотя и тоже очень легко читаемый, чем «бед» файл. «Бед» файл – крайне прост. У нас есть хромосома, которая может быть не обязательно хромосома, это может быть «контиг сф», вот, которая просто именует фрагмент, скажем, мы начинаем дальше, через «таб» у нас должны быть две еще обязательно колонки, то есть, а потом есть количество необязательных колонок. Но, опять же, нам это пока неважно, вот. И у нас есть начало, то есть, и конец, то есть, у нас должны быть координаты начала и конца в «бед» файле, точно также можно записать не только гены, можно записать снипы, указав, ну то есть, у нас начало, конец совпадают любые «инделы». На самом деле, любую информацию в геноме можно описать интервалами. То есть, если вы хотите более изощренного какого-то описания, то есть, ну, на самом деле, все, что угодно. То есть, «бед» файл можно усложнять и в том числе, описывать там интронные зоны и все остальное. Вот, но в самом простом варианте, он действительно крайне примитивен и его можно сделать хоть руками. Но, мы этого делать не будем. А мы попробуем сейчас посмотреть, насколько файлы, которые у нас получилось сделать использовать «глиммер» и «джинмарк» подходят под описание «бед» файла. Итак, еще раз проверяем, что у нас здесь. Смотрите, то есть у нас есть bacteria_genome.fasta.lst, давайте мы его скопируем на уровень выше, то есть... Оп! И смотрим на уровень выше. Он у нас тут есть? Да, он у нас тут есть. И давайте сюда теперь скопируем из глиммера аннотацию, которую с координатами. Отлично! Скопировали. У нас здесь как-то все замусорилось, давайте мы немножко отсортируем то, что получилось. Создадим временную директорию TMP, в которую сложим, на самом деле нам индексы уже точно сейчас не понадобятся, архивные .vcf-файлы мы тоже уберем, осталось только glimmer_germs.coords, и sorted_bacteria.bam давайте тоже прячем, вот. И продолжаем работать. Проверяем, что у нас все на месте. И как же нам вообще все это сделать? На самом деле, используйте любой текстовый редактор. Я предлагаю использовать TextWrangler, он бесплатный и доступный. И начнем с того, что у нас сделано – bacteria_genome.fasta. Ага. Смотрим, что получилось. В общем-то, как мы видим, у нас здесь есть вся информация, которая нам нужна. У нас есть, то есть хромосома здесь у нас одна и та же, вот она, она у нас в заголовке стоит, а дальше она нигде не указывается, но мы это можем легко добавить в строки. Дальше, что у нас есть? У нас есть тренд, который, в общем-то, информация полезная, но не критичная сейчас, она необязательна в .bed-файле, но давайте мы ее добавим, у нас есть же необязательные поля. У нас есть левый и правый хвост файла. Отлично. Размер гена и все остальное уже не так обязательно. Как же мы это можем сделать? На самом деле, еще одну вещь давайте посмотрим. Мы посмотрим, какими, то есть у нас в .bed-файле они должны быть разделены табами, а здесь у нас чем разделены между собой строки? Для этого надо показать невидимые символы, о! На самом деле, это большой кошмар! Они разделены пробелами, причем не одним пробелом, а большим количеством пробелов. Это крайне неприятно, но поправимо. Можно, в общем-то, все делать руками, то есть взять и все строки отредактировать, их там 3000, но мы так делать не будем. Давайте мы попробуем сделать иначе. Мы сейчас воспользуемся регулярными выражениями в простом виде. Они встроены, на самом деле, во все системы, наверное, какие только есть, и использовать их тоже крайне просто, вы зачастую их используете и даже не осознаете, то есть всякие варианты поиска и так далее, вот. И в TextWrangler'е надо открыть поиск и отметить, что мы используем grep. Отлично! Как же он работает? Давайте посмотрим! То есть вот если мы ищем просто символ «1», мы легко это делаем, мы можем искать плюс. Ага, только плюс, поскольку он в регулярных выражениях имеет определенное значение, мы просто плюс не можем найти, смотрите: мы можем искать как бы символ «+». Отлично, что мы можем делать еще? На самом деле, а дальше начинаются интересные варианты. Одним из простых, а мы попробуем найти... мы можем создавать свои условия поиска, то есть в квадратных скобках мы можем указать то, что мы можем «или, или», то есть мы будем удовлетворены, если это будет плюс или минус, смотрите, и оно отлично работает. Нам это, возможно, где-нибудь потом пригодится, вот. И нам нужно каким-то образом заменить разное число пробелов, потому что вы видите, что у нас значения разной длины и пробелов между строками разное на одинаковые табы. Как мы это можем сделать? Смотрите, у нас есть замечательный символ начала строки, который определяет, что у нас, ну, с чего начинается строка. А дальше есть «\s» – означает пробел, вот мы всегда находим по одному пробелу с начала строки. Но мы можем это сделать все еще лучше: мы можем сказать «\s+», и это будет означать– найти максимальное число пробелов после начала строки. Команда, в общем-то, крайне нетрудная и легко выполнимая. И опять же в принципе все это, все редактирование можно было б сделать с помощью различных методов, я просто вам хочу показать, как это можно сделать с помощью grep'а, причем я буду стараться использовать крайне простые выражения, чтобы вам было понятно, как оно работает. И мы нажимаем теперь, смотрите: менять это надо просто, мы убираем символы пробела с начала строки, нажимаем replace all. Отлично, смотрите! У нас уже изменился файл. У нас пропали пробелы в начале, то есть это мы сделали, теперь. Теперь нам нужно каким-то образом убрать вот эти вот интервалы из пробелов между колонками. Воспользуемся «\s», ставим плюс, смотрите-ка, да, пробелы у нас находятся. Давайте теперь попробуем сказать, что мы интервалы из пробелов заменяем «\t», это обозначение таба, то есть символ, обозначающий таб. Меняем. О! Все чудовищно испортилось. Меняемся назад. В чем же беда? Давайте попробуем поискать. Ага, оказывается, смотрите, такое выражение, в том числе, считает, что у нас конец строки тоже содержит серию пробелов, и это неприятно. Как же нам это сделать? На самом деле, вариантов-то, наверно, много, но мы с вами сейчас страшно схитрим. Смотрите, мы его заменим на какой-то символ, который у нас в файле не встречается, например, плюс-минус; и заменим не все такие, а мы целенаправленно изменяем только конец строки. В конце строки в линуксовых системах у нас, просто конец строки мы заменим на другой символ, да, вот. Это «\n», означающий конец строки. В Windows'е он должен быть там «\r», «\n», вот. Но, я надеюсь, никто Windows'ом сейчас не пользуется, потому что это вам будет чересчур трудно, если вы попытаетесь повторить. Вот. И мы меняем файл. У нас все испортилось, кажется, но на самом деле мы теперь можем сказать, что у нас, мы меняем все, все пробелы, серии пробелов мы меняем на табы. Замена, а после этого мы восстанавливаем то, что мы сделали: мы опять несуществующие символы в файле меняем на конец строки. Ух-ты! Смотрите, великолепно! Мы уже переформатировали файл, но на этом у нас не заканчивается вся мощь регулярных выражений, потому что .bed-файл-то мы так и не сделали. Мы просто привели файл в форму, которая нам достаточно удобна. Так, что же нам нужно будет сделать. Нам теперь нужно сформировать будет такие строки, в которых будет хромосома, а она у нас на самом деле всегда одна, вот она, то есть нам нужно каким-то образом будет вставлять ее туда, и координаты начала и конца. Начало и конец – это у нас, получается, третья колонка и четвертая, вот. И сейчас я покажу другую вещь, связанную с регулярными выражениями, которая тоже является крайне мощной, смотрите. Мы, используя «\d», которое ищет любой цифровой символ, можем, собственно, легко это сделать. «\d+» означает, что мы ищем самый длинный. То есть после этого мы ставим, что мы ищем «\t», как видите, он у нас находит уже длиннее, после «\t» у нас в нашей строке будет всегда идти два символа, точнее либо один, либо другой: это будет либо минус, либо плюс, за которым опять идет «\t». Отлично, как видите, мы можем все время проверять, что у нас он правильно находит, но у нас нашлась не вся строчка, мы дальше ищем «\d», опять «\t» и опять «\d+». И, да, это выглядит не столь красиво, когда мы хотим это всю строку сделать, да. Но в конце у нас, кстати, смотрите, в конце поскольку у нас всегда будут идти единицы, то мы можем опять же поставить сюда просто единицу, смотрим. Мы теперь несмотря на то, что эти все строки у нас разные, вот это паттерн, который мы задали в строке поиска, он совпадает со всеми из перечисленных. Ну, по-моему, это уже достаточно удивительно, но если вам этого мало, а вам этого должно быть мало, потому что .bed-файл мы так и не сделали, то мы делаем следующую вещь. Смотрите, мы хотим, то есть выделяя круглыми скобками квадратные скобки, у нас это, мы формируем свой как бы набор символов для поиска, который будет. Круглые скобки – это что-то вроде запоминания. И можно, если мы скажем, что у нас круглая скобка идет от, начиная от «\t» до символа плюс-минус, дальше у нас идет от «\t» до первого «\d+», точнее, это будет третья колонка, и четвертая колонка, вот. Таким образом, мы с вами, как бы, сейчас выделили круглыми скобками колонку номер два, три, четыре. Отлично! А теперь мы воспользуемся командой replace в новом варианте, смотрите. На самом деле, когда мы делаем вот это выделение круглыми скобками, мы как бы запоминаем переменную, и мы ее теперь можем использовать, потому что у нас по порядку они идут, то есть у нас есть три переменных, которые так и запоминаются сейчас: это будет «\1», «\2» и «\3». И нам нужно сначала «\2», потому что мы плюс добавим, то есть у нас начало, «\3» и «\1». Смотрим, что у нас получилось. Удивительное дело, смотрите! Мы сформировали самый замечательный красивый файл с tab delimited, то есть он у нас соответствует спецификации .bed, у нас есть хромосома, начало и конец, и ориентация на стренде. Отлично! На самом деле, вся остальная часть нам уже здесь даже и не нужна. Прекрасно получилось. Собственно, сделали мы это весьма быстро, то есть вам не пришлось писать никаких сложных программ, действительно, и это делается все просто. И на самом деле, это очень полезно, потому что очень часто приходится работать с переходами из одного формата в другой и так далее. То есть регулярные выражения – это, можно сказать, то, из чего можно, по крайне мере, в некоторых способах использования, это то, с чего стоит начать учить, вот. Нам осталось только сохранить файл. Давайте мы его сохраним как, это будет, мы же это считаем, что у нас .bed-файл, и так его и назовем, то есть это у нас будет .bed. И, сейчас, это у нас он исходный был, это был genmark, поэтому это будет genm_arks.bed, genm_marks.bed, .bed-файл. Отлично, сохраняем. Ну что, у нас здесь появилось genm_marks.bed? Проверяем. Ой. Отлично. Вот, с этим разобрались.
