Что же делать дальше? Дальше нам нужно разобраться со вторым файлом, который у нас сделал glimmer. Открываем glimmer. Может там будет всё лучше? Нет, здесь, кажется, то же самое, что у нас было и тогда. Здесь у нас тоже опять всё в пробелах, которые нам совершенно не годятся, и к тому же ещё и нет заголовка, который описывает. Но, во-первых, полей в общем-то меньше и вообще всего меньше, то есть в чём-то лучше, но, на самом деле, работать будем практически так же. Спрятали невидимые буквы, они нам не нужны. Нам потребуется, нам потребуется имя хромосомы, здесь оно не указано, и давайте опять воспользуемся регулярными выражениями. Смотрите, у нас есть все орфы начинаются со слова, собственно, orf и дальше набор цифр. Мы его так и можем искать, смотрите. И мы все орфы вполне благополучно давайте заменим на. . на хромосому. Бац, сделали. Прекрасно. Символов пробелов в начале строки нет, так что нас это никак не пугает. Давайте прибегнем к тому же страшному костылю, который делался до этого у нас. Есть /n, потому что я подозреваю, что так он тоже не сработает. Меняем его на. Оп! Дальше мы меняем серию пробелов на табы и возвращаем символы конца строки на место. Угу, отлично, всё получилось. На самом деле можно, можно здесь использовать регулярные выражения, но зачастую это может быть даже трудно и, может, неоправданно. И я вам покажу, что можно использовать ещё. То есть у нас поскольку здесь колонки даже сделаны в общем-то в том виде, в котором они нам нужны, то мы можем сказать, что это просто glimmer точка, – давайте назовем его glimmer_pre.bed, то есть это у нас не совсем .bed файл, но мы где-то на пути к нему. Смотрим, есть ли у нас эти все файлы тут. Да, glimmer_pre.bed, вот он. Давайте заглянем ещё раз внутрь файла, что у нас корректно открывается. Прекрасно. У нас есть в терминале замечательная команда cut, которая у нас умеет резать, и смотрите как это можно сделать. Мы говорим, мы f минус один и указываем тот же самый glimmer_pre.bed, только использовать, – файл у нас большой, поэтому если вы не хотите, чтобы у вас всё вылетело в терминал, – используйте команду head: он покажет тогда только макушку. Смотрите, мы получили просто первую колонку. Если мы скажем, что мы хотим, – здесь, кстати, пробелы необязательно, – мы можем сказать: мы хотим первую, вторую и третью. Опа! Прекрасно, получили то, что хотели. Вот, дальше мы точно так же можем это сохранить в новый файл, использовав символ перенаправления. Если вы хотите, кстати, сейчас... Просто не всегда разделяются табами. Если бы у вас разделение шло пробелами или вообще любыми другими символами, я, кстати, могу даже продемонстрировать. Смотрите. Мы ставим, мы хотим первую колонку, но при этом мы укажем, что у нас делимитер – это точка. Смотрите, он, по сути, ищет первую точку в файле и по ней делит, как бы, всё. Но нам нужно не это, а нам нужно было, чтобы вот так. И теперь записываем в новый файл glimmer.bed. Вот. Ну как мы видим, здесь у нас получился он чуть покороче, чем в предыдущий раз. То есть у нас нет ориентации, но, вроде как, и бог с ней. Сохраняем. Будьте очень аккуратны, работая в терминале, потому что на самом деле вот символы записи файла или удаления, ещё что-то, то есть у вас здесь редко спрашивают, уверены ли что вы хотите это сделать, и вы можете несколькими неловкими командами убить дни своей работы или часы. В общем, будьте аккуратны и внимательны, когда что-то делаете. Заглядываем внутрь. Получилось. Проверим, что у нас получилось. По размеру он тоже никак не изменился, нет, всё прекрасно. Ну а теперь всё просто: мы сделаем пересечение этого файла, то есть intersectBed. Работает. У нас должен быть первый файл glimmer и у нас должен быть второй файл genemarks. Угу, смотрите! Он не очень сработал, и говорит, что у нас неправильный .bed файл. И действительно, мы же были не очень внимательны. На самом деле я, в общем-то, специально оставил, потому что просто показать те ошибки, которые вы можете встретить. В .bed файле у нас должны быть отсортированы начало и конец, а у нас получается в этой строке: третья координата, ну то есть вторая координата, третья колонка меньше, чем первая координата во второй колонке. И это надо как-то менять. Давайте посмотрим на файл. Посмотрим на файл, который glimmer.bed. Ну и да, вот, например, это такая строка. Видите, у нас начало гена, а потом после. Что же делать? А вот тут опять же вариантов есть самое большое множество, но я хочу вам показать, как эту проблему можно решить с помощью крайне небольшого скрипта на «Питоне». Я знаю, что программирование на «Питоне» не входит в текущий курс, но я надеюсь, что, во-первых, простота того, что там происходит, и то, как это можно использовать, сподвигнут многих из вас на то, чтобы начать учить Python, потому что без него всё равно далеко в биоинформатике не продвинуться, вот. И именно с помощью «Питона» мы сейчас попробуем выяснить, а где же ориентация у этих генов и каким-то образом их рассортировать. И для демонстрации возможностей «Питона» я приготовил маленький скрипт, который попытается решить ту проблему, которая у нас возникла. Собственно, под названием «Где же моя ориентация», в текущем случае – «наша». Смотрите. Как видите, здесь всего есть 20 строк, которые решают нашу проблему. Из того, как бы, что здесь происходит, алгоритм, который вам подсказывают, здесь как раз вот его и видно. Что мы берём каждую строку файла, – вот здесь говорится for line in в файле, который мы ему дадим, – мы ломаем эту строку сначала, то есть мы должны со строки убрать символы конца строки, которые /n мы искали грепом, они нам всё равно не нужны, и дальше ломаем его на куски. Куски, кусков у нас получается три, и мы эти три куска записываем в переменную outList, которая будет просто списком из этих трёх кусков. А дальше мы должны сравнить координаты. Если у нас, то есть опять же в «Питоне» и вообще в программировании у нас получается, часто индексация идёт, по-моему, чаще всего с нуля, а не с единицы. То есть для тех, кто вообще этим не занимался, поначалу будет надо привыкнуть, то есть первый символ в списке – он имеет порядковый номер ноль, а потом уже идущие за ним будет один и два. То есть у нас есть три элемента: ноль, один и два. Вот, мы сравниваем элемент первый с координаты один с элементом с координаты два. Если он меньше, то есть у нас получается начало координаты меньше координаты конца, – всё правильно. Значит все, что нам нужно сделать, это добавить к нашему списку символ плюса и конца строки, чтобы у нас файл корректно писался потом, и мы на этом останавливаемся. Если это не так, то есть в случае, когда у нас не выполняется вот это условие, у нас происходит другая вещь – значит у нас, наоборот, значит у нас сначала идёт конец, потом начало. И в этом случае с помощью вот такой вот команды мы можем поменять местами, то есть мы просто говорим две переменных и указываем, что дальше они меняются местами. И пришиваем к хвосту в этом случае опять же минус, то есть мы говорим, что у нас ориентация будет обратная. После всего этого мы формируем новую строку с помощью из списка outList и команды join, используя здесь у нас стоит «таб», мы формируем новую строку, которую после этого пишем в файл. Файл после этого закрываем. Как видите, здесь всё крайне то есть инструктивно, я имею ввиду, что вам, как вы думаете, что надо делать, как вы бы это, может, руками делали, так тут и происходит, но просто не руками, а всё сделает «Питон». И да, собственно, вот это будет glimmerСorrected.bed, по идее это будет имя файла, который получится. И если вы будете писать сами свой файл, то перед тем, как его запустить, вам нужно будет дать разрешение на его запуск. Сделать это командой chmod с такими параметрами, а дальше давайте попробуем посмотреть, что у нас получится. И мы ему даем glimmer, glimmer.bed Вот. У нас отработала корректно, по крайней мере, никаких ошибок не выдала. Давайте смотрим, что внутри. У нас есть glimmerCorrected.bed Отлично. Смотрим, как он выглядит внутри. Выглядит внутри он очень симпатично. Смотрите. У нас, во-первых, всегда теперь, — по крайней мере тут, где мы видим, — начало меньше конца, плюс добавился символ с ориентацией. То есть, в общем-то, его и так было понятно, но теперь мы можем на самом деле делать различные совершенно вещи, например, мы можем, используя grep -c и поиск по символу «минус», мы можем в этом файле найти все гены и пересчитать их, которые имеют обратную ориентацию, чего мы до этого не могли. И таких приблизительно половина. Собственно, без использования параметра «минус» мы можем просто получить строки этих файлов, если вы, например, хотите заниматься именно ими по какой-то причине. Вот. Отлично. У нас получилось вроде как два корректных bed-файла, которые мы теперь, наконец-то, можем использовать с intersectBed. Давайте, кстати, посмотрим параметры интерсекта еще раз. У него есть важный параметр -f, который говорит о том, какой у нас оверлэп должен быть между двумя интервалами. По дефолту это единица, то есть перекрытие на один нуклеотид. Для нас, конечно, это очень мало. То есть, если бы мы делали может другие вещи типа снипов, было б нормально, но гены у нас должны совпадать лучше. Давайте они у нас будут совпадать на 90 %. После этого мы говорим, что у нас glimmer, glimmerCorrected.bed и genemarks.bed. Отлично. Они пересекаются. То есть мы видим, что ответ получается и, да, сколько у нас всего таких получилось. Мы это уже сделаем с помощью полюбившейся команды word count и строки. Смотрите-ка. У нас на 90 % совпали результаты glimmer и genemarks. То есть не на 90 %, а с перекрытием 90 % у нас совпали предсказания, что, в общем-то много. То есть крайне большое количество генов не совпало. Мы можем уменьшить строгость, например, до 80 %, и количество будет увеличиваться. То есть при 70 оно будет еще выше. Как вы видите, обе программы выдали достаточно хорошие и сходящиеся результаты, что на самом деле очень хорошо, но будьте готовы к тому, что на других биологических объектах, где модель, возможно, не столь хорошо подходит, у вас могут быть значительные отклонения в том, что находят разные тулы. И это было всего лишь две программы, а программ по поиску генов великое множество, и каждая из них может в разных условиях дать что-то свое. Вот. Но у нас получился очень неплохой результат, в котором у нас совпали, совпали не только предсказания у программ между собой, они совпали по своему количеству с тем, что проаннотировано в оригинале, то есть 3000... около 3000 белок-кодирующих генов. Вот. То есть, в общем-то, у нас есть аннотация генов, дальше мы можем вернуться к нашим... к нашим снипам, и посмотреть снипы попадают ли у нас в гены, и, вообще, что там происходит. Для этого есть опять же специальные инструменты типа SNP-эффекта, который можно установить, сделать к нему базу или скачать готовую, если она есть. И он сделает все за вас. Или вы можете воспользоваться тем же самым intersect.bed-ом. Параметры может vcf-файл воспринимать. И посчитаем просто, в какие, в какие гены у нас попадает сколько, сколько снипов. Вот. В этом случае нам нужно взять будет... -b мы меняем на var_calls.vcf. Отлично. -f нам уже не нужно. У нас снипы — это одна буква, но нам нужен другой параметр. И, если я не ошибаюсь, это -c. Да. Он будет считать, сколько же у нас раз пересеклось с каждым... с каждым интервалом. Для каждого интервала он сделает дополнительную колонку, в которой покажет, смотрите, сколько у нас пересеклось с ней генов. Он показывает 0. И как вы думаете, почему у нас ничего не получилось, что тоже бывает предсказуемо. Смотрим на наш vcf-файл. На самом деле 98 % всех ваших трудностей будут связаны с форматами. И смотрите, мы же использовали просто имя хромосомы, когда... в bed-файле для генов, а здесь у нас имя хромосомы совсем другое. Оно включает не только там... Оно и NCBI-ный айдишник и genebank-овский. Вообще, оно длиннее. И программа его узнать не может. И в этом случае вы правы совершенно, если думаете о том, что мы воспользуемся мощью регулярных выражений для того, чтобы исправить эту плачевную ситуацию. И мы должны будем искать фрагмент, который начинается с gi. Плюс — это любой символ. То есть мы ищем любые символы до тех пор, — максимально длинную строку, — пока не встречаем вертикальную черту, на которой у нас заканчивается это. И... О! Мы таким образом паттерным выхватили полностью имя в vcf-файле скаффолда. Ну, и меняем просто его на то, к которому мы привыкли. Отлично. У нас удалось отредактировать vcf-файл. Теперь сохраняем и возвращаемся к тому, что мы делали. И смотрим еще раз. О! Теперь программа может узнать интервалы между собой и все нам посчитала. И как вы видите, мы отсортировали сразу даже еще и по количеству снипов в каждом из генов. Таким образом, у нас... удивительно, но у нас есть 3 гена, содержащие огромное количество по сравнению с другими... то есть под сотню снипов, то есть в отличие... Если у всех предыдущих их них там десятки, а у большинства будут нули. Вот. Как вы видите, можно получить очень много информации из того, что получилось, потому что вот эти гены можно смотреть вручную. По-хорошему, дальше нужно продолжать аннотировать получившиеся варианты, используя поиск по... разделив их на категории синонимичность и не синонимичность и так далее, но это делается уже сложнее, и я думаю, что... То есть опять же вы можете это сделать, если сделаете специальную базу для SNP-эффекта или напишите полностью свой скрипт. Но в текущем варианте, я думаю, это уже более чем достаточно, и мы с вами научились многому. Мы мало того, что нашли гены. Мы смогли скомбинировать эту информацию с разных источников и мы пересекли ее с теми данными, которые получили на прошлой неделе и нашли три самых вариабельных участка, предсказанных как генами в геноме исследуемой бактерии. Спасибо за внимание.
