Этап 1

Для работы использовали бактерию Vescimonas fastidiosa . Белковые последовательности его протеома были загружены из раздела UniProt Proteomes, а через эту запись получены идентификаторы соответствующей геномной сборки.

ID протеома в UniProt Proteomes: UP000681035
Код доступа геномной сборки из GenBank: GCA_018408575.1
Код доступа геномной сборки из RefSeq: GCF_018408575.1

Этап 2

Чтобы скачать последовательность генома и таблицу локальных особеностей использовалась данная команда:
datasets download genome accession GCF_018408575.1 --include genome,gff3

--include - уточняет, какие конкретно файлы в составе указанного генома нужны
genome - запрос геномной последовательности
gff3 - запрос таблицы локальных особенностей

Архив с файлами был распакован командой: unzip ncbi_dataset.zip

Этап 3

Определение варианта генетического кода организма, выполненное с помощью загрузки данных таксона из NCBI Taxonomy, стало первым этапом перед поиском открытых рамок считывания:
efetch -db 'taxonomy' -id '2714355' -format 'xml'

efetch - извлекает записи из указанных баз данных в виде отдельных файлов
-db - опция указания базы данных
-id - опция указания TaxID
-format - эта опция определяет формат вывода данных
2714355 - Taxid V. corpula в NCBI

В поле был указан №11 генетического кода, этот номер часто встречается у бактерий, архей и пластид растений

Поиск открытых рамок считывания и их трансляций: getorf -sequence ncbi_dataset/data/GCF_018408575.1/GCF_018408575.1_ASM1840857v1_genomic.fna -outseq orf.faa -minsize 150 -table 11 -find 0

-sequence указание на файла на вход с последовательностью
-outseq - имя выходного файла
-find 0 - поиск рамок между стоп-кодонами
-minsize 150 - минимальное значение количества нуклеотидов, принимаемое на вход как открытую рамку считывания.
-table 11 - указание номера генетического кода

Нужно проверить, что среди трансляций нет тех, которые короче 50 а.о.:
infoseq orf.faa -only -length | sort -n -u | less
-only -length - выводит только длины последовательностей
sort -n -u - сортировка по возрастанию длины и вывод уникальных значений
less - чтение полученного

Создадим белковую базу для blastp командой:
makeblastdb -in orf.faa -dbtype prot -out proteome
-in - входной файл
-dbtype - указанине типа последовательности
-out - выходной фай

Этап 4

Получить последовательности гомологичных метилтрансфераз возможно командой:
echo "sw:P0AED9 sw:P0AEE8 sw:P23941" | tr ' ' '\n' | seqret @stdin -osformat fasta -outseq query.fasta
@stdin - читает входные данные из стандартного ввода

Этап 5

Для выравнивания трёх последовательностей с транслированными ORF из генома использовали команду:
blastp -query query.fasta -db proteome -out outcome -outfmt 7
-outfmt 7 - формат выдачи (табличная)

Выдача blastp

Наибольшая по весу находка имеет индификатор NZ_AP023418.1_4753 (123)
Координаты: [950508 - 951416] (команда: cat orf.faa | grep 'NZ_AP023418.1_4753' )
Гомологична m6A-МТазе E.coli

Определение СDS из таблицы локальных особенностей генома, которые располагаются рядом. Сначало получим СDS из той же геномной последовательности:
grep -E '^NZ_AP023418.1.*CDS' ncbi_dataset/data/GCF_018408575.1/genomic.gff | cut -f4,5,7,9 > CDS.tsv
cut -f4,5,7,9 - извлекает столбцы
Теперь отберем из этого файла СDS c близками координатами к находке:
echo -e '950508\t951416\t-\tFOUND-ORF' | cat - CDS.tsv | sort -n | grep -C 3 'FOUND-ORF' > neighbors.tsv
echo -e '950508\t951416\t-\tFOUND-ORF' - в поток передается строка: 950508‹TAB›951416‹TAB›-‹TAB›FOUND-ORF
cat - CDS.tsv - принимает ввод из предыдущей команды (строку с FOUND-ORF) и объединяет эту строку с содержимым файла CDS.tsv
grep -C 3 'FOUND-ORF' - показывает 3 строки до и 3 строки после совпадения и ищет строку, содержащую FOUND-ORF

Выдача

В файле была обнаружена CDS (WP_213541931.1), координаты которой (950520 – 951419) частично перекрывают координаты нашей находки (950508 - 951416). В продуктах также была указана ДНК-метилтрансфераза (product=Dam family site-specific DNA-(adenine-N6)-methyltransferase) что является прямым подтверждением гомологии исследуемой последовательности.

Этап 6

Проведем поиск по аннотациям кодирующих участков:
elink -db 'nuccore' -id 'NZ_AP023418.1' -target 'protein' | efilter -query '2.1.1.72[EC/RN Number]' | efetch -format 'gp'
-db nuccore - указывает исходную базу данных Nucleotide
elink - для поиска связанных записей между разными базами данных NCBI
efilter -query '2.1.1.72/2.1.1.32/2.1.1.113[EC/RN Number]' - на вход поступает список всех id белков из шага 1, а на выходе остается только множество этих id, которые соответствуют ферментам (имеют EC number).
gp - GenPept полная запись о белке

По m6A конвейер выдал 1 запись: WP_213541931.1
По m4C и m4C ничего не нашел записей.