Геном Aeromonas allosaccharophila

Изучаем геном Aeromonas allosaccharophila и ищем ДНК-метилтрансферазу.

AC

Чтобы получить AC геномной сборки, надо сначала получить TaxID организма, которому принадлежит геном:

zcat aer_allo.swiss.gz | grep '^OX' | cut -b 1-20 | sort | uniq -c | cat

4008 OX   NCBI_TaxID=656

Получим список сборок, соответствующих таксону; я решил воспользоваться конвейером с esearch и efetch:

esearch -db 'taxonomy' -query ''656' [UID]' | elink -target 'assembly' | efetch -format 'uid'

Полученная выдача - список номеров сборок.

Получим список AC по данному таксону:

datasets summary genome taxon 656 --as-json-lines | dataformat tsv genome | cut -f1 | sort | uniq | less

Полученная выдача - список AC. GCA_016026615.1 - нужный. Убедимся, что он есть в списке:

grep 'GCA_016026615.1' 656.ac

GCA_016026615.1

Убедились, всё в порядке.

Загрузка генома и таблицы локальных особенностей

Загрузим нужные таблицы с помощью полученного ранее AC:

datasets download genome accession GCA_016026615.1 --include gff3 --include genome

Распакуем:

unzip ncbi_dataset.zip

Поиск ORF

Для начала надо узнать, какая таблица генетического кода у моей бактерии:

efetch -db taxonomy -id '656' -format xml | less

Полученная выдача - xml таблица

grep -A 1 '<GeneticCode>' 656.tax # 656.tax - файл, в который я сохранил выдачу

    <GeneticCode><GCId>11</GCId>

По выдаче видно, что нужна таблица 11.

Получим ORF:

getorf -sequence GCA_016026615.1_ASM1602661v1_genomic.fna -table 11 -minsize 150

Промежуточный файл

Сделаем базу данных blast:

makeblastdb -in 656.orf -out 656.bla -dbtype prot

Проверим их длину:

infoseq -sequence 656.orf -only -length | sort -n | uniq | less

Убедился, длина более 50 нуклеотид.

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

У ДНК-метилтрансфераз прокариот есть ряд гомологических доменов. Воспользуемся для поиска ДНК-метилтрансферазы Aeromonal allosaccharophila последовательностями белков-гомологов - m5C-MTase и m6A-MTase E. coli и m4C-MTase Bacillus amyloliquefaciens:

echo sw:P0AED9 sw:P0AEE8 sw:P23941 | tr ' ' '\n' > s.txt | seqret @s.txt query_MTases.fasta

По сделанной ранее БД ORF применим blastp:

blastp -task blastp -query query_MTases.fasta -db ORFs -out blastp.out -evalue 0.05 -outfmt 7

Выдача - таблица XML в формате 7

Получено всего две находки со значимым e-value - CP065745.1_33961 (6.71e-122) по запросу DMA_ECOLI P0AEE8 DNA adenine methylase и CP065746.1_12 (7.29e-16) по запросу MTB1_BACAM P23941 Type II methyltransferase M.BamHI. Для DCM_ECOLI P0AED9 DNA-cytosine methyltransferase ничего не было найдено.

Отдельно перезапишем координаты и дополнительную информацию CDS нашей геномной последовательности:

grep 'CDS' genomic.gff | grep 'CP065745.1' | cut -f 4,5,7,9 > CDS.tsv

Выдача

Посмотрим, какие CDS пересекаются с моей находкой.

echo -e '2368072\t2368995\t-\tMY-ORF' | cat - CDS.tsv | sort -n | grep -C 3 'MY-ORF' > neighbors.tsv

Полученные соседи

Моя находка целиком покрывает QPR52938.1. Любоптыно, что первый триплет CDS не попадает в мою рамку считывания (старт CDS - 2368069, старт моей рамки считывания - 2368072), и что ORF заканчивается на 54 нуклеотида дальше от конца CDS (концы: CDS - 2368941, ORF - 2368995). Продукт - Dam family site-specific DNA-(adenine-N6)-methyltransferase

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

Попробуем найти CDS находки по EC-коду фермента:

elink -target protein -db nuccore -id CP065745.1 | efilter -query '2.1.1.72' | efetch -format acc > out.txt

Выдача состоит из одного AC - QPR56285.1. По какой-то причине QPR52938.1 не был найден.