EMBOSS, Entrez Direct, NCBI Datasets

Этап 1: получение АС геномной сборки

Получение TaxID выбранного мной в прошлом году организма - Cupriavidus cauae:

zcat UP000324324.swiss.gz | grep 'TaxID=' |  cut -c 17-23 | sort -u | less

Выдача команды: 2608999

Получение списка сборок, соответствующих таксону:

datasets summary genome taxon 2608999 --as-json-lines | dataformat tsv genome | vd

В выдаче 2 сборки: одна для моего организма Cupriavidus cauae - GCF_026210475.1, а другая для организма Cupriavidus gilardii - GCF_008632125.1. Я, соответственно, выбрала первую.

Этап 2: скачивание файлов

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

datasets download genome accession GCF_026210475.1 --include gff3,genome

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

Определение типа генетического кода:

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

Ответ: тип 11 (стандартный для бактерий)

Нахождение рамок считывания и получение их трансляций (длина которых при этом не менее 50 аминокислот):

getorf -filter 'GCF_026210475.1_ASM2621047v1_genomic.fna' -table 11 -minsize 150 getorf.fasta

Создание локальной белковой базы для blastp с названием ORFs:

makeblastdb -in getorf.fasta -dbtype prot -out ORFs

Проверка длины трансляций:

infoseq -only -length getorf.fasta | sort -u -n | less

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

Скачивание последовательности белка ДНК-метилтрансфераза прокариот: P0AED9 (Dcm, m5C-МТаза, E.coli), P0AEE8 (Dam, m6A-МТаза, E.coli), и P23941 (m4C-МТаза, Bacillus amyloliquefaciens) 

echo 'sw:P0AED9' 'sw:P0AEE8' 'sw:P23941'| tr ' ' '\n' > query_1MTases.fasta 

seqret @query_1MTases.fasta query_MTases.fasta

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

Я провела blastp по созданной базе ORFs в поисках сходств последовательностей ДНК-метилтрансфераз у упомянутых выше бактерий:

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

Файл с табличной выдачей.

Нашлось два сходства. Лучшая по весу находка - первая, для m5C-МТаза E.coli.

Таблица 1. Выдача локального blastp
№ рамки у моей бактерии Координаты рамки у моей бактерии Bit score % identity Обнаруженный гомолог
3195 652471 - 653772 364 51.122 Dcm

По координатам находки blastp я нашла какие CDS из таблицы локальных особенностей генома располагаются рядом:

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

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

Файл с выдачей. Получилось, что найденная мной открытая рамка пересекается с двумя CDS, но координаты не совпадают: моя рамка имеет координаты 652471 - 653772, а CDS dcm (на + цепи) 652507 - 653775 и CDS (на - цепи) 652074 - 653105

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

Я решила проверить, можно ли найти CDS, полученные в предыдущем пункте, с помощью поиска по аннотации кодирующих участков в геноме:

elink -target 'protein' -db 'nuccore' -id 'NZ_CP080293.1' | efilter -query '2.1.1.37' | efetch -format 'acc'

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