Учебная страница курса биоинформатики,
год поступления 2023
Задания практикума 9
Предполагается, что все программы вы будете запускать на kodomo, не требуется установка каких-либо программ на свой компьютер.
При желании, некоторые из программ вы можете установить себе.
EMBOSS доступен для Linux, Windows (хоть и не самая последняя версия), и, кажется, Mac. Инструкции и ссылки для скачивания доступны тут. На kodomo есть установочный файл для Windows, если у вас не получается справиться с FTP-сервером sourceforge.
EDirect – на сайте есть только версия для Linux и Mac.
Datasets CLI – существуют версии для всех ОС, есть даже github.
Blast+ доступен для всех OC, инструкции здесь.
Проверять буду страницу с отчетом на сайте kodomo. Она может быть в виде либо HTML, либо Wiki, как вам удобнее. Ссылку на страницу нужно обязательно указать при записи на проверку, и где-нибудь на своем сайте. Мягкий дедлайн – 01:00 AM 5 ноября, жесткий дедлайн – 01:00 AM 12 ноября. После мягкого дедлайна штраф 0.5 балла, после жесткого – 2 балла.
Настоятельно рекомендую первым делом создать папку ~/term3/pr9, и все делать внутри нее (можно в подпапках). Не в public_html! В public_html скопируйте только файлы, которые требуются в отчете, чтобы привести на них ссылки со страницы. После сдачи практикума папку ~/term3/pr9 можно будет удалить, освободив место, но не испортив отчет. Промежуточных файлов может быть много, но проверять их я не буду. Не замусоривайте домашнюю папку!
Глобальная задача: с помощью средст EMBOSS, EDirect, NCBI Datasets CLI и blast+ найти одну из ДНК-метилтрансфераз в геноме бактерии/археи по последовательности и по аннотации. Заодно придется освежить в памяти UniProt и bash.
Имейте в виду, в первую очередь от вас требуется освоить технические средства. Поэтому в отчете должны быть приведены все использованные команды. Код должен быть оформлен подобающим образом (как минимум, шрифт должен быть моноширинным). Однако текст отчета должен быть понятен в отрыве от заданий. Не пропускайте вводную информацию и пояснения. Выдачу команд приводить в отчете не нужно, кроме случаев, когда это явно требуется в задании.
Не буду проверять страницы, которые не удовлетворяют перечисленным выше формальным требованиям. За каждую перепроверку из-за несоответствия формальным требованиям буду снимать 0.5 балла. Из перепроверок по другим причинам штрафуется только первая.
Копировать любой текст из заданий или, тем более, со страниц других студентов строго запрещается. Обращаться за советом или пояснениями можно и нужно, но все команды и все пояснения вы должны сформулировать самостоятельно.
Этап 0: протеом бактерии/археи
Вы будете исследовать тот самый прокариотический геном, на основе которого получен протеом, который вы сравнивали с контролем в прошлом семестре (практикум 8). Вам потребуются белки этого протеома, скаченные в формате swiss. Они могли сохраниться с прошлого семестра, в этом случае этот этап у вас уже выполнен. Если вы удалили файл для экономии места, то найдите в отчете команду, которую использовали для скачивания, и повторите её.
Этап 1: получение AC геномной сборки
Можно было бы получить эту информацию с сайта UniProt, но этот способ сложнее автоматизировать. Можете использовать его для проверки, в отчете это упоминать не нужно.
Для этого сначала нужно получить TaxID организма, которому принадлежит геном. В формате swiss для этой информации предусмотрено соответтсвующее поле. C помощью grep и других консольных программ обработки текста (tr, cut, sort, ...) получите содержимое этого поля для всех белков протеома и убедитесь, что все значения совпадают. Приведите конвейер и выдачу в отчете.
В базе NCBI Assembly (которая в веб-интерфейсе теперь называется Datasets Genome) многие сборки присутствуют в двух версиях: GenBank (GCA_...) и RefSeq (GCF_...). Версии RefSeq, как правило, лучше аннотированы, чем GenBank (в крайнем случае, это просто полная копия). Используйте версию RefSeq в случае её наличия.
Получить список сборок, соответствующих таксону, можно как с помощью программы datasets, так и с помощью конвейера esearch + efetch. Используйте любой вариант. Результатом команды, приведенной в отчете, должен быть список сборок в виде таблицы, содержащей, кроме прочего, столбец с AC. Не надо сразу скачивать последовательности или таблицы локальных особенностей.
Если в списке несколько сборок (не являющихся разными версиями одной сборки), то вам придется определить, какая из них соответствует вашему протеому. Для этого найдите в белках протеома ссылки на нуклеотидные записи INSDC (вспоминайте практикум 7 из прошлого семестра). По списку этих AC получите соответствующие записи в Assembly с помощью elink. Это должно сократить список до одной сборки, если вы все сделали правильно. Приведите все команды и пояснения в отчете.
Этап 2: скачивание последовательности генома и таблицы локальных особенностей
Имея AC геномной сборки, последовательность и feature table легко загрузить с помощью команды datasets download genome accession, добавив нужные аргументы и опции. Правда таблицу локальных особенностей можно загрузить только в формате GFF3, он значительно отличается от того, что вы использовали в первом семестре. Но информации в нем содержится даже больше. Читайте справку, полное описание формата приведено здесь. Но то, что вам понадобится, вполне можно найти и без полного понимания формата.
После загрузки архив с файлами надо распаковать с помощью программы unzip.
Этап 3: поиск и трансляция открытых рамок считывания
Прежде, чем искать открытые рамки считывания, нужно определить, какой вариант генетического кода использует исследуемый организм. Эту информацию можно получить из записи про таксон в базе NCBI Тaxonomy, скачав её с помощью efetch в формате xml. Информация про таблицу кода содержится и в базе Nucleotide, можно скачать summary записей про геномные последовательности, входящих в сборку (с помощью esummary). Выбирайте любой вариант. Не все бактерии и археи используют таблицу №11, может быть вам повезло выбрать организм с каким-нибудь альтернативным генетическим кодом.
Найти открытые рамки считывания и сразу получить их трансляции можно с помощью программы getorf из пакета EMBOSS. Вам потребуется получить рамки между стоп-кодонами (чтобы не потерять фрагмент последовательности из-за неверного старт-кодона), при этом их трансляции должны быть не короче 50 аминокислот, чтобы совсем маловероятные короткие фрагменты отбросить (их будет очень много). Не забудьте указать нужную таблицу кода! Последовательности трансляций сохраните в какой-нибудь промежуточный файл. Создайте по ним белковую базу для blastp, назовите её ORFs.
С помощью программы infoseq из EMBOSS вам надо обязательно проверить, что среди трансляций нет тех, которые короче 50 а.о. Вы могли неверно указать порог длины рамок считывания. В отчете надо привести проверочную команду без выдачи. Выдачу просмотрите сами.
Вместо предварительного поиска рамок и трансляции можно было бы использовать tblastn. У обоих подходов есть плюсы и минусы, сами подумайте, в чем они состоят.
Если в будущей работе вам потребуются не только трансляции, но и нуклеотидные последовательности рамок считывания, то можно попросить getorf выдать нуклеотиды (опция -find), а трансляции получить с помощью программы transeq.
Этап 4: получение последовательностей гомологичных метилтрансфераз
Все известные ДНК-метилтрансферазы прокариот, по всей видимости, содержат гомологичные каталитические домены. Однако они могут быть насколько далеки друг от друга, что сходство последовательностей стандартными средствами обнаружить не удается. Тем не менее, скорее всего у вас получится найти как минимум одну ДНК-метилтрансферазу по сходству последовательностей со следующими белками (указаны коды доступа в базе Swiss-Prot): P0AED9 (Dcm, m5C-МТаза, E.coli), P0AEE8 (Dam, m6A-МТаза, E.coli), и P23941 (m4C-МТаза, Bacillus amyloliquefaciens).
Скачать последовательности этих белков из Swiss-Prot можно с помощью seqret. Постарайтесь сделать это одним конвейером, вспомнив про специальную форму USA, которая позволяет использовать сразу несколько адресов, записанных в файл. Содержимое для этого файла легко генерируется командой конвейером из echo и tr, причем можно это сразу передать в seqret без создания промежуточного файла, надо только придумать, как заставить seqret читать listfile из STDIN.
Файл с последовательностями назовите query_MTases.fasta
Этап 5: поиск по сходству последовательностей
С помощью локального blastp, база ORFs, query_MTases.fasta в качестве запроса. Ссылку на табличную выдачу в формате 7 надо привести в отчете. Дальше анализируйте только лучшую по весу находку. В отчете нужно привести название рамки, её координаты в геноме, вес находки, гомолог какой из МТаз был обнаружен (m5C, m6A или m4C). Если у вас есть сомнения по поводу гомологичности находки, то их надо описать в отчете.
По координатам находки blastp можно определить, какие CDS из таблицы локальных особенностей генома располагаются рядом. Для этого сначала надо отобрать нужные строки и столбцы из таблицы локальных особенностей и записать их, например, в файл CDS.tsv. Строки нужны те, которые соответствуют CDS в той же геномной последовательности, что и ваша находка, а столбцы – с координатами (4 и 5), цепью (7) и дополнительной информацией (9). Обратите внимание, что в таблице локальных особенностей в столбце 4 всегда меньшая координата, не всегда это начало кодирующего участка. Этим можно воспользоваться, а именно добавить строчку, соответствующую вашей находке, отсортировать записи по столбцу меньшей координате, а потом отобрать с помощью grep соседние CDS.
echo -e '12345\t14567\t-\tMY-ORF' | cat - CDS.tsv | sort -n | grep -C 3 'MY-ORF' > neighbors.tsv
Ссылку на neighbors.tsv надо привести в отчете. По содержимому этого файла вам надо определить, пересекаются ли какие-нибудь CDS с найденной вами открытой рамкой. Если да, то совпадают ли координаты.
Этап 6: поиск по аннотациям кодирующих участков
В завершении, вам надо определить, смогли бы вы найти CDS, соответствующий (-ие) вашей находке, с помощью поиска по аннотации кодирующих участков в геноме.
Поиск по аннотации нужно производить с помощью конвейера из elink (по АС нуклеотидной записи в nuccore получаете белковые записи в protein), efilter (отбираете только те белки, которые имеют подходящие аннотации), efetch (получаете AC отобранных белковых записей).
Поиска нужно провести по EC-кодам ферментов (для ДНК-метилтрансфераз это 2.1.1.37 (m5C), 2.1.1.72 (m6A) и 2.1.1.113 (m4C)). Выберите тот, который соответствует вашей находке, или ищите сразу по трем. В отчете кроме конвейера edirect надо привести количество найденных белков и указать, присутствует ли среди них тот, кодирующий участок для которого в аннотации генома пересекается с найденной рамкой считывания.
Использовать логические операторы в запросах efilter можно только вместе с группирующими скобочками, иначе он может сломаться и начать выводить записи далеко не только из фильтруемого набора. Возможно, исходный набор добавляется к запросу в виде дополнительного условия через AND без проверки, что запрос остается правильным.
Поиск по названию белка в базе NCBI Protein производить сложно, так как можно искать только целые слова. Например, если белок называется "Site-specific DNA-methyltransferase", то по слову "methyltransferase" вы его не найдете. Для поиска по названиям белков используйте UniProt.