EMBOSS, Entrez Direct, NCBI Datasets

Этап 1. Получение AC геномной сборки

Для получения AC необходимо знать TaxID организма. Его можно получить при помощи конвеера. Необходимая информация содержится в файле после в строке, начинающейся с OX, каретка перед названием строки необходима, чтобы отсечь все строки, где сочетание букв OX встречается внутри текста.

zgrep ^OX UP000054387.swiss | tr -s ' ' | cut -f2 -d ' ' | sort -u

Вывод: NCBI_TaxID=1514971

AC был получен при помощи команд NCBI Datasets:

datasets summary genome taxon 1514971 --as-json-lines | dataformat tsv genome

Для моего таксона нашлась только одна геномная сборка.
AC: GCF_001469955.1

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

Feature table и последовательность можно скачать с помощью команды из NCBI Datasets:

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

Формат GFF3 представляет собою модифицированный GFF. Файлы GFF3 состоят из девяти столбцов, разделенных табуляцией. Один объект в GFF3 (например экзон) может принадлежать более чем к одной группе одновременно, имеется строгое разграничение имени/идентификатора элемента и предположения о его принадлежности к той или иной категории.

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

Вначале необходимо определить вариант генетического кода, используемый Haloprofundus marisrubri.

efetch -id '1514971' -db 'taxonomy' -format 'xml' > code.xml

Из xml видно, что это стандартный вариант — таблица номер 11. Следовательно, в getorf указываем опцию -table 11. Чтобы ограничить длину рамки считывания снизу, используем опцию minsize. Команда:

getorf GCF_001469955.1_ASM146995v1_genomic.fna file.fasta -minsize 150 -table 11 -find 0 -filter

Убедимся, что всё сработало верно и в файле действительно отсутствуют последовательности меньше 50 а.к.:

infoseq file.fasta -length -only | sort -hr

Создаём новую базу данных ORFs:

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

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

Скачивание последовательностей метилтрансфераз из Swiss-Prot производится с помощью конвеера:

echo 'sw:P0AED9' 'sw:P0AEE8' 'sw:P23941'| tr ' ' '\n' | seqret -filter '@stdin' query_MTases.fasta

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

Осуществим поиск:

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

Находки: файл
Наибольший вес имеет находка NZ_LOPU01000029.1_4790 (гомолог m5C-МТазы): 63.2, E-Value: 4.63e-11. Координаты: [831863 - 833062]
Далее работаем с ней.
Необходимо определить CDS из таблицы локальных особенностей генома, которые расположены рядом. Отберём нужные строки и столбцы из таблицы и запишем их в файл CDS.tsv.

grep 'NZ_LOPU01000029.1' genomic.gff | grep 'CDS' | cut -f4,5,7,9 > CDS.tsv

Определим участки, координаты которых лежат рядом с координатами находки:

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

Выдача файла доступна здесь
Просмотрев выдачу, можно обнаружить, что находке лучше всего соответствует CDS с координатами [831872 - 833065] (почти полностью совпадают). Из аннотации следует, что это ген DNA (cytosine-5-)-methyltransferase, то есть мы обнаружили гомолог m5C-МТазы.

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

Осталось попробовать найти CDS, соответствующий моей находке, с помощью поиска по аннотации кодирующих участков в геноме. Сделаем это с помощью конввера, а затем оценим результат:

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

Программа выдала 1 результат. Он соответствует белку, кодирующий участок которого в аннотации генома пересекается с найденной рамкой считывания.

Источники информации:

Haloprofundus marisrubri proteome