EMBOSS, Entrez Direct, NCBI Datasets

Петренко Павел

Факультет биоинженерии и биоинформатики, Московский Государственный Университет имени М.В.Ломоносова

1. Протеом бактерии

В 8-ом практикуме прошлого семестра я исследовал протеом бактерии Moraxella bovoculi 237 с идентификатором UP000035860. Он у меня уже был скачан, поэтому я его просто скопировал в директорию term3/pr9.

2. Получение AC геномной сборки и TaxID организма

3. Скачивание последовательности генома и таблицы локальных особенностей

Теперь, имея АС геномной сборки, я скачиваю последовательности генома и таблицы локальных особенностей с помощью команды:

datasets download genome accession GCF_000696305.2 --include gff3,genome #скачал архив
unzip ncbi_dataset.zip #распаковал архив

В результате, в директории ncbi_dataset/data я получил два json-файла: dataset_catalog.json с информацией о всех полученных файлах и assembly_data_report.jsonl с информацией сборки в NCBI (как я понял), а также папку GCF_000696305.2 с необходимыми мне файлами: GCF_000696305.2_MBO_237_genomic.fna (последовательность генома в формате fasta) и genomic.gff (таблица локальных особенностей в формате GFF3).

4. Поиск и трансляция открытых рамок считывания

Сначала на NCBI Taxonomy Browser я взял ID для Moraxella bovoculi 237: 743974. Затем я определил какой вариант генетического кода использует Moraxella bovoculi 237 с помощью команды efetch (лично мне efetch кажется более простой в использовании, к тому же она даёт больше данных):

efetch -db taxonomy -id '743974' -format 'xml'> taxonomy.xml #скачал информацию и записал в файл taxonomy.xml

В файле нашли нужную строку и выяснили, что бактерия использует таблицу №11:

<GeneticCode>
<GCId>11</GCId>
<GCName>Bacterial, Archaeal and Plant Plastid</GCName>
</GeneticCode>

p.s. Могли просто посмотреть номер таблицы генетического кода на сайте NCBI Taxonomy Browser.

С помощью getorf из EMBOSS находим открытые рамки считывания:

getorf ncbi_dataset/data/GCF_000696305.2/GCF_000696305.2_MBO_237_genomic.fna translation.fna -minsize 150 -table 11 -circular -find 0

Несколько пояснений про код: -circular делает так, что ORF может искаться в участках соединения начала и конца последовательности (то есть он как бы замыкает последовательность), -find 0 транслируем области между стоп-кодонами (что и требуется по заданию).

Теперь infoseq из EMBOSS проверим, что все транляции не короче 50 аминокислот:

infoseq translation.fna -only -length | sort -n > aa_length #-only -length выводит только длины последовательностей

Посмотрев файл aa_length с длинами транслированных последовательностей, убедились, что они не короче 50 аминокислотных остатков. Создадим по полученным трансляциям белковую базу для blastp и назовём её proteome:

makeblastdb -in translation.fna -dbtype prot -out proteome

5. Получение последовательностей гомологичных метилтрансфераз

На этом этапе надо скачать указанные в задании ДНК-метилтрансферазы прокариот, для этого воспользуемся конвеером из команд echo, tr и seqret:

echo "sw:P0AED9 sw:P0AEE8 sw:P23941" | tr ' ' '\n' | seqret -filter 'list::stdin' -outseq query.fasta

Несколько пояснений про код: echo "sw:P0AED9 sw:P0AEE8 sw:P23941" - записываем в строку нужные нам коды доступа, tr ' ' '\n' - переносим по строкам (меняем пробелы на переносы строк), -filter 'list::stdin' - подаём файл из stdin (-filter вообще используем всегда, хуже не будет, как нам сказали на паре).

6. Поиск по сходству последовательностей

С помощью blastp был осуществлён поиск гомологов указанных раньше метилтрансфераз по созданной белковой базе данных бактерии Moraxella bovoculi 237:

blastp -query query.fasta -db proteome -outfmt 7 -out blastp_result

Файл с табличной выдачей результата работы blastp (формат 7)

Всего получилось 13 находок (4 m5C, 2 m6A и 7 m4C). Далее рассмотрим гомолог m6A, являющийся лучшим по весу среди всех находок.

Теперь найдём, какие CDS располагаются рядом. Для этого сначала отберём нужные нам строки из таблицы локальных особенностей с помощью конвеера:

grep 'NZ_AOMT01000023.1' ./ncbi_dataset/data/GCF_000696305.2/genomic.gff | grep 'CDS' | cut -f 4,5,7,9 > CDS.tsv

Несколько пояснений про код: grep 'NZ_AOMT01000023.1' - отобрали строки с нужными AC, grep 'CDS' - отобрали строки с нужными кодирующими последовательностями, cut -f 4,5,7,9 - оставили нужные столбцы (координаты, +/- цепи и дополнительная информация).

Поищем соседей для лучшего выравнивания blastp. Сначала найдём координаты последовательности в файле с открытыми рамками считывания, затем с помощью конвеера, указанного в задании, найдём "соседей":

grep 'NZ_AOMT01000023.1_274' translation.fna #нашли координаты

Видим, что координаты (92734 - 93618) расположены по возрастанию, следовательно у нас + цепь.

echo -e '92734\t93618\t+\tFOUND-ORF' | cat - CDS.tsv | sort -n | grep -C 3 'FOUND-ORF' > neighbors.tsv

Файл с выдачей

Найденно перекрывание с геном WP_036365547.1, который расположен на + цепи - это тоже m6A-МТаза (как и наша находка). Также есть ещё одно частичное перекрывание с некоторым предсказанным белком на - цепи (WP_036365550.1), перекрывание составляет лишь 5 нуклеотидов, что может возникнуть из-за оперенной организации или регуляторных перекрываний.

7. Поиск по аннотациям кодирующих участков

Поищем CDS моей сборки (NZ_AOMT01000023.1) по аннотациям. Для этого сделаем конвеер из предложенных в задании команд (elink, efilter и efetch), сравнивать будем с 2.1.1.72 (m6A), соответствующей моей находке:

elink -target 'protein' -db 'nuccore' -id 'NZ_AOMT01000023.1' | efilter -query '2.1.1.72' | efetch -format 'fasta'

Несколько пояснений про код: elink -target 'protein' -db 'nuccore' -id 'NZ_AOMT01000023.1' - нашли белковые фрагменты в базе protein, транслируемые сборками, похожими на исходную, efilter -query 2.1.1.72 - фильтруем результат по нужной метилтрансферазе, efetch -format 'fasta' - получили белковые последовательности в формате fasta, автоматически сохранились в файле query.fasta.

В результате видим тот самый белок, который мы нашли в предыдущем пункте - WP_036365547.1, который соответствует m6A-МТазе.