Практикум 9

1 Получение AC геномной сборки

В первую очередь надо получить AC геномной сборки, а для этого я должен получить TaxID моей бактерии из файла UP000182840.swiss.gz
Для этого я выполню следующий конвейер:

Выдача: NCBI_TaxID=1670800
Пояснение: в результате работы конвейера остался только один уникальный TaxID, значит у всех белков в файле они одинаковые.

Затем я получил список геномных сборок при помощи следующей команды:
datasets summary genome taxon 1670800 --as-json-lines | dataformat tsv genome
Либо так, чтобы удобнее было смотреть:
datasets summary genome taxon 1670800 --as-json-lines | dataformat tsv genome | vd
Сборок несколько, и они разные (у них разные Assembly name), так что надо выбрать ту, которая соответствует имеющемуся у меня протеому. Я посмотрел на штамм в файле протеома и в информации о сборках и пришёл к выводу, что надо выбрать GCF_001889605.1, так как там тот же штамм, что и в протеоме (B7), а ещё это версия RefSeq, то есть она лучше аннотирована.

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

Теперь скачиваю последовательность генома и таблицу локальных особенностей, используя полученный AC и программу datasets. Для этого выполняю команду:
datasets download genome accession GCF_001889605.1 --include gff3,genome --filename hahagenomefiles.zip
Далее распаковываю архив: unzip hahagenomefiles.zip

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

Прежде чем искать открытые рамки считывания, нужно узнать, какой генетический код использует изучаемый организм. Эту информацию можно найти в базе данных NCBI Taxonomy с помощью следуюущей команды:
efetch -db 'taxonomy' -id '1670800' -format 'xml' | less
Из строк, содержащихся в выдаче (Рис. 1.), следует, что Aquibium oceanicum использует стандартную для бактерий таблицу генкода № 11.

Рис. 1. Информация о таблице генкода

Следующим моим шагом стал поиск открытых рамок считывания в геноме и их трансляция. При этом продукты трансляции должны быть не короче 50 аминокислот. Такие короткие я не учитываю, чтобы отмести огромное количество маловероятных продуктов.
Команда:
getorf GCF_001889605.1_ASM188960v1_genomic.fna translated.fasta -table 11 -minsize 150
Пояснения: 150 — минимальная длина открытой рамки считывания, значит минимальная длина пептида будет 150/3=50.
GCF_001889605.1_ASM188960v1_genomic.fna — входной файл с последовательностью генома, а translated.fasta — выходной файл.

Затем я решил проверить точно ли не получилось последовательностей меньше 50 нуклеотидов.
Команда: infoseq -only -length translated.fasta |sort -n -u
Пояснения: использую infoseq чтобы вывести только длины транслированных аминокислотных последовательностей, а затем сортирую их по возрастанию и оставляю только уникальные значения для удобства.
Теперь можно создать белковую базу для blastp: makeblastdb -in translated.fasta -out ORFs -dbtype prot

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

Известно, что прокариотические метилтрансферазы вероятно содержат гомологичные каталитические домены, но сильно различаются по последовательности. Я попробую найти в Aquibium oceanicum метилтрансферазу, похожую на одну из следующих: P0AED9 (Dcm, m5C-МТаза, E.coli), P0AEE8 (Dam, m6A-МТаза, E.coli), и P23941 (m4C-МТаза, Bacillus amyloliquefaciens). Используя USA форму и коды доступа в базе Swiss-Prot я скачал нужные последовательности:
1. echo -e 'sw:P0AED9\nsw:P0AEE8\nsw:P23941' > mt.txt
2. seqret @mt.txt query_MTases.fasta

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

Проведу blastp по базе ORFs и query_MTases.fasta в качестве запроса:
blastp -task blastp -query query_MTases.fasta -db ORFs -out blastp_hehe.out -outfmt 7 Результат в табличной форме
По результатам поиска нашлось 6 результатов для m5C-, 3 для m6A- и 6 для m4C-метилтрансферазы.
Лучшая находка была обнаружена для m4C-метилтрансферазы — NZ_CP018171.1_7704, ее вес(bit_score) = 89.0, а координаты: 1487256 — 1488437 (координаты нашёл в файле translated.fasta)
Далее следовало определить, какие кодирующие участки (CDS) пересекаются с найденным мною участком. Я создал файл CDS.tsv, содержащий CDS только из той же нуклеотидной последовательности, что и найденный участок, а также дополнительную информацию:
grep 'NZ_CP018171.1' genomic.gff | grep 'CDS' | cut -f4,5,7,9 > CDS.tsv
Затем были получены участки близкие к целевой области:
echo -e '1487256\t1488437\t-\tMY-ORF' | cat - CDS.tsv | sort -n | grep -C 3 'MY-ORF' > neighbors.tsv
neighbors.tsv
Да, получается что как минимум одна CDS(1487307-1488440) пересекается с моим участком(1487256-1487256), но координаты не совпадают. Но из аннотации следует, что эта CDS - DNA-methyltransferase, то есть ДНК-метилтрансфераза, то есть как раз гомолог m4C-метилтрансферазы. Класс!

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

Напоследок я поставил перед собой задачу проверить можно ли найти эту CDS, пересекающуюся с найденным мной участком при помощи поиска по аннотации кодирующих участков в геноме:
elink -target protein -db nuccore -id NZ_CP018171.1 | efilter -query '2.1.1.113' | efetch -format 'fasta'
Пояснение: искал я используя EC-код нужной мне метилтрансферазы (m4C).
Увы, найти таким образом ничего не удалось.

7 СОПРОВОДИТЕЛЬНЫЕ МАТЕРИАЛЫ

1. Справка по getorf