Ну и да, вот, например, это такая строка.
Видите, у нас начало гена, а потом после.
Что же делать?
А вот тут опять же вариантов есть самое большое множество, но я хочу вам показать,
как эту проблему можно решить с помощью крайне небольшого скрипта на «Питоне».
Я знаю, что программирование на «Питоне» не входит в текущий курс,
но я надеюсь, что, во-первых, простота того, что там происходит, и то, как
это можно использовать, сподвигнут многих из вас на то, чтобы начать учить 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-файл.
