Для работы с элементами генома организма необходимо скачать геномную сборку, а для этого нужно в первую очередь определить TaxIDганизма. Для этого к файлу протеома был применен следующий конвейер
grep "^OX" UP000001013.swiss | cut -d ' ' -f 4 | tr -d ';'| sort -u
NCBI_TaxID=186497
Далее по уже известному идиентификатору были найдены относящиеся к данному организму геномные сборки. С помощью следующего конвейера их было найдено 2 различных сборки, каждая в 2 версиях - GenBank и RefSeq.
datasets summary genome taxon 186497 --as-json-lines | dataformat tsv genome
Для того, чтобы понять, которая из них относится к данному протеому найдем по идеинтификаторам нуклеотидных записей, к которым относятся белки, в базе нуклеотидных последовательностей ссылки на геномную сборку:
grep "DR" UP000001013.swiss | grep 'EMBL;' | cut -d ' ' -f 5 | tr -d ';'| elink -db 'nuccore' -target 'assembly' | esummary | xtract -pattern DocumentSummary -element AssemblyAccession
В выдаче была одна геномная сборка: GCF_000007305.1
Скачаем геномную сборку и таблицу локальных особенностей к ней и разархивируем.
datasets download genome accession GCF_000007305.1 --include genome --include gff3
unzip ncbi_dataset.zip
Перед тем как искать ORF необходимо определить таблицу генетического кода данной архею, на случай если она вдруг не является стандартной. Для этого по ID организма скачаем данные о некоторых параметрах организма.
efetch -db 'taxonomy' -id '186497' -format xml
Выснено что код организма описывается таблицей 11. Теперь мы можем при помощи пакета EMBOSS найти возможные ORF. Тк код триплетный, а минимальные ORF мы полагаем не меньшими 50 аминокислотных остатков, ставим минимум в 150 нуклеотидов.
getorf -sequence 'GCF_000007305.1_ASM730v1_genomic.fna' -table 11 -minsize 150 -out porf
Проверка на ошибку: выделим только длину полученных последовательностей и отсортируем по возрастанию и возьмем верхушку. Наименьшие среди полученных последовательностей - 50 а.о.
infoseq -only -length porf | sort -V -u| head
Создадим базу из полученных ORF средствами локального BLAST, чтобы после по ней найти гомологи к ферментам среди возможных последовательностей.
makeblastdb -in porf -dbtype prot -out ORFs
Скачаем последовательности метилаз E.coli, чтобы потом найти к ним гомологов.
echo sw:P0AED9,sw:P0AEE8,sw:P23941 | tr ',' '\n' | seqret -filter '@stdin' query_MTases.fasta
Для поиска гомологов применил алгоритм blastp скачанным последовательностям. Наибольший вес был у последовательности у рамки считывания NC_003413.1_11395 и составлял 59.3, гомология для нее определилась с m5C-МТазой. Координаты этой ORF в геноме 56557 - 57459
blastp -query query_MTases.fasta -db ORFs -outfmt 7 -out MTase_blast
Определим по таблице локальных особенностей, пересекается ли эта последовательность с каким-нибудь CDS. Для этого из таблицы возьмем строки относящиеся к CDS и к тому же участку в геноме, отсортируем их по координатам, выбрав для выдачи нужные нам столбцы - координаты и некоторые особенности.
grep 'NC_003413.1' genomic.gff | grep 'CDS' | cut -f4,5,7,9 > CDS.tsv
echo -e '56557\t57459\t-\tMY-ORF' | cat - CDS.tsv | sort -n | grep -C 3 'MY-ORF'
Выяснено, что наша последовательность пересекается с двумя CDS. С одной (55848-56570) только самым краем, но со второй (56554-57459) совпадает практически полностью, только координаты начала сдвинуты на 3. При этом в свойствах есть запись:DNA (cytosine-5-)-methyltransferase activity, что свидетельствует о том, что определенная последовательность действительна гомологична. Идиентификатор - WP_011011161.1.
Посмотрим, можно ли прийти к такому же результату через поиск по аннотации кодирующих участков. Искать будем по тому же участку генома. Проводим поиск по EC коду фермента - m5C метил-трансферазы.
elink -db nuccore -target protein -id 'NC_003413.1' | efilter -query '2.1.1.37' | efetch -format fasta
В выдаче ровно одна последовательность белка - WP_011011161.1, это то же самый, что нашли до этого - DNA cytosine methyltransferase.