В 8 практикуме прошлого семестра мною был выбран и скачан протеом UP000000799. Основная задача данного практикума состоит в том, чтобы найти с помощью технических средств EMBOSS, EDirect, Datasets CLI ген(-ы) ДНК-метилтрансфераз в геномной сборке, соответствующей этому протеому.

• Страница сборки на UniProt Proteomes — Campylobacter jejuni subsp.
• Идентификатор протеома — UP000000799
• Код доступа геномной сборки из GenBank - GCA_000009085.1
• Идентификатор RefSeq - GCF_000009085.1
• NCBI Taxonomy ID - 192222

a) Имея AC сборки, скачиваем нуклеотидную последовательность генома и таблицу локальных особенностей с помощью программы:

• datasets download genome accession GCF_000009085.1 --include gff3 --include genome

--include gff3(genome) — получить файл с последовательностью генома и таблицу особенностей (feature table) в формате gff3

b) Скачанный архив ncbi_datasets.zip был разархивирован следующей командой:

• unzip ncbi_dataset.zip

a) Прежде чем транслировать полученную нуклеотидную последовательность нужно было выяснить какую таблицу генетического кода использует мой организм, поиск был проведен следующей командой:

• efetch -db 'taxonomy' -id '192222' -format 'xml' > tax.xml
• -db — база данных, в которых проводится поиск
• -id — TaxID бактерии
• -format — формат выходного файла (xml)

В полученном файле tax.xml было обнаружено, что данный организм использует таблицу №11. Информация о используемой бактерией генетической таблице находится в поле GeneticCode.

b) Поиск открытых рамок считывания был осуществлен с помощью следующей программы:

• getorf GCF_000009085.1_ASM908v1_genomic.fna trans.fasta -minsize 151 -table 11 -find 0
• GCF_000830845.1_ASM83084v1_genomic.fna – файл со скачанной нуклеотидной последовательностью генома
• trans.fasta – файл с искомыми открытыми рамками считывания
• -minsize 151 – минимальная длина ORF в нуклеотидах
• -table 11 – используем таблицу генетического кода №11
• -find 0 – искать ORF только между стоп-кодонами, игнорируя старт-кодон (чтобы не потерять фрагменты из-за неизвестного старт-кодона).

* в задании указано фильтровать ORF меньше 50 аминокислот. каждая аминокислота кодируется 3 нуклеотидами → 50 * 3 = 150 нуклеотидов. 151 означает чуть больше 50 а.о., чтобы включить границу

c) Далее была осуществлена проверка, что среди транслированных белков нет тех, что короче 50 а. о., с помощью программы:

• infoseq trans.fasta | tail -n +2 | cut -f3 -d '-' | sort -n -k2,2

Вывод программы infoseq показал, что все полученные трансляции являются белковыми (тип P) и имеют длину не менее 233 аминокислот.

• infoseq trans.fasta – выводит таблицу с информацией о трансляциях
• tail -n +2 – пропускаем заголовок, берём только данные
• cut -f3 -d '-' – берём третий элемент, разделённый дефисом
• sort -n -k2,2 – сортируем по числовому значению второго поля, чтобы увидеть диапазоны длин ORF

d) Далее была создана белковая база данных на основе полученных открытых рамок считывания:

• makeblastdb -in trans.fasta -dbtype 'prot'
• -dbtype 'prot' – указываем, что база белковая

В результате выполнения команды makeblastdb была создана белковая база данных BLAST, представленная набором служебных файлов (.psq, .pin, .phr и др.)

Теперь необходимо было скачать последовательности гомологичных метилтрансфераз (по сходству с которыми можно было бы осуществить поиск генов этих белков у данного организма), а именно:

1. Dcm, m5C-МТаза, E.coli (код доступа в Swiss-Prot P0AED9).
2. Dam, m6A-МТаза, E.coli (код доступа в Swiss-Prot P0AEE8).
3. m4C-МТаза, Bacillus amyloliquefaciens (код доступа в Swiss-Prot P23941)

Скачивание последовательностей осуществлялось в 2 этапа, сначала был создан listfile tr.txt:

• echo -e 'sw:P0AED9\nsw:P0AEE8\nsw:P23941' > tr.txt

Далее была применена программа seqret и создан файл query.fasta, содержащие нужные последовательности:

• seqret @tr.txt query.fasta

a) Наконец, с помощью blastp, созданной ранее базы данных, а также скачанных в прошлом пункте последовательностей в качестве запроса (fasta-файл query_MTases.fasta), был произведен поиск гена метилтрансферазы у данного организма:

• blastp -task blastp -query query.fasta -db trans.fasta -out blast.out -outfmt 7
• blastp - программа из пакета BLAST+, которая сравнивает белковую последовательность (protein) с белковой базой данных
• -task blastp - тип алгоритма
• -query query.fasta - файл с запросными белковыми последовательностями
• -db trans.fasta - имя BLAST-базы данных (использует все файлы на выдаче)
• -out blast.out - файл, в который будет записан результат поиска
• -outfmt 7 - 7 табличный формат с комментариями

Ознакомиться с выдачей blastp можно здесь. Самая лучшая находка по весу, как можно видеть из выдачи, имеет идентификатор NC_002163.1_972, самый высокий bit score (31.2), самый маленький e-value (0.11), координаты в геноме 581824 - 584013.

команда для поиска координат:

• grep '^>NC_002163.1_972' trans.fasta
b) Далее нужно было определить с какими кодирующими белки участками (CDS) пересекается моя лучшая находка, для этого сперва нужно было создать вспомогательный файл CDS.tsv, содержащий только CDS из той же нуклеотидной последовательности, что и моя находка, а также необходимую информацию о них (главным образом координаты и аннотация), все данные о CDS были взяты из ранее скачанной таблицы локальный особенностей (файл genomic.gff)
• grep NC_002163.1 genomic.gff | grep CDS | cut -f4,5,7,9 > CDS.tsv

Далее были обнаружены участки, которые по координатам наиболее близки (или даже перекрываются) с моей находкой:

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

В файле neighbors.tsv лежит выдача, с которой можно ознакомиться здесь. Отсюда видно, что найденная находка очень хорошо соответствует CDS с координатами 581827 - 584016. Индефикатор YP_002344052.1(carbamoyltransferase). Однако это не метилтрансфераза

В завершении был опробован также поиск нужных CDS по аннотациям в геноме, для чего были использованы EC-коды соответствующих ферментов: 2.1.1.37 (m5C), 2.1.1.72 (m6A) и 2.1.1.113 (m4C). Поиск осуществлялся конвейером:

• elink -target protein -db nuccore -id NC_002163.1 | efilter -query '2.1.1.113' | efetch -format 'fasta'
• elink ищет связи между записями в разных базах данных NCBI.
• -target protein — указывает, что мы хотим найти связанные белки (Protein)
• -db nuccore — указывает исходную базу данных (нуклеотидные последовательности)
• efilter -query '2.1.1.113' - фильтрует результаты поиска
• efetch -format 'fasta' - получает данные из базы NCBI в fasta-формате

К сожалению, ничего не нашлось. Это объясняется особенностями работы метода. elink + efilter ищет только официально аннотированные белки, связанные с данной нуклеотидной записью, и фильтрует их по EC-номеру. Если белок присутствует в геноме, но не имеет официальной аннотации EC в базе NCBI, он не будет найден этим способом.