EMBOSS, Entrez Direct, NCBI Datasets

Протеом бактерии

Для практикума 8 в прошлом семестре мы работали с референсным протеомом Bartonella bovis 91-4 , который включает 1267 последовательностей белков. Скачаем его с помощью команды: 

wget 'https://rest.uniprot.org/uniprotkb/stream?compressed=true&format=txt&query=%28%28proteome%3AUP000014038%29%29' -O UP000014038.swiss.gz

Получение AC геномной сборки и TaxID организма

Скачали необходимую информацию используя АС для RefSeq (последовательность и feature table):

datasets download genome accession GCF_000384965.1 --include genome --include gff3 --filename 91-4.zip

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

Со страницы NCBI Genomes мы взяли NCBI Taxonomy ID — 1094491 и использовали команду efetch чтобы найти тип генетического кода:

efetch -db taxonomy -id '1094491' -format 'xml'> taxonomy.xml

В выдаче нашли соответствующую строку:

<GeneticCode>
  <GCId>11</GCId>
</GeneticCode>

Затем мы транслировали открытые рамки считывания с помощью getorf. 

 getorf -sequence ./ncbi_dataset/data/GCF_000384965.1/GCF_000384965.1_BBbMcIrM01_genomic.fna -minsize 150 /
 -table 11 -circular -find  -outseq open_frames.fasta
Нам нужны результаты трансляции не менее 50 аминокислот, поэтому установили параметр -minsize 150 (длина параметра указывается в нуклеотидах)
-table 11 — указание генетического кода
-circular — анализ кольцевой последовательности
-find — трансляция рамок между стоп-кодонами (по умолчанию это параметр 0)

Проверим, что правильно отсортировали результаты трансляции по длине >50 аминокислот. Для этого выведем их по возрастанию длины с помощью команды из пакета EMBOSS: 

 infoseq -sequence open_frames.fasta  -only -length|sort -n|uniq|less
Первая строка в выдаче — 50, то есть меньше нет.

Построим по найденным открытым рамкам белковую базу для blastp: 

makeblastdb   -in open_frames.fasta   -dbtype prot   -title proteome   -out proteome
— программа создала 8486 последовательности

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

Мы хотим среди полученных открытых рамок считывания попробовать угадать белок метилтранферазу. Для этого протестируем наши последовательности на гомологичность другим метилтрансферазам. Для загрузки последовательнотей трёх претендентов: P0AED9 (Dcm, m5C-МТаза, E.coli), P0AEE8 (Dam, m6A-МТаза, E.coli), и P23941 (m4C-МТаза, Bacillus amyloliquefaciens) был составлен конвейер.
Три AC для swissprot подаются как список, затем его можно подать seqret двумя эквивалентными способами (@mylist и list::mylist). Ну и фильтр добавим чисто для кайфа. 

echo 'sw:P0AED9 sw:P0AEE8 sw:P23941' | tr ' ' '\n' | seqret -filter 'list::stdin' -outseq query.fasta

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

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

Мы сделали поиск по blastp — на вход даём три белковые последовательности метилтрансфераз и сравниваем с нашей приготовленной белковой базой данных. 

blastp -query query.fasta -db proteome -out MTase.out -outfmt 7
Выдача конвейера с выравниваниями

Всего нашлось 23 хита для всех трёх метилтрансфераз. У всех высокий E-value, кроме одной находки среди хитов для MTB1_BACAM P23941 и её m4C-МТазы:

Поищем последовательность, к которой выровнялась метилтрансфераза среди открытых рамок считывания, найденных ранее: 

 grep '^>NZ_CM001844.1_1078' open_frames.fasta
— это фрагмент генома [464413 - 465696]. Заметим, что координата начала меньше координаты конца, то есть фрагмент на (+) цепи.

Соседние кодирующие последовательности к найденной гомологичной

Для поиска соседних CDS сначала отберём их от всех последовательностей. Отбираем те, которые начинаются с нужного AC и только CDS, затем нужные колонки: 4,5 - координаты, 7- направление цепи, 9 - дополнительная информация. Записываем в файл. 

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

Затем мы берём найденные координаты лучшего выравнивания blastp и выводим 3 CDS перед ней и 3 после, чтобы посмотреть на соседей.

echo -e '464413\t465696\t+\tFOUND-ORF' | cat - CDS.tsv | sort -n | grep -C 3'FOUND-ORF' > neighbors.tsv

Выдача конвейера с соседними CDS

Часть выдачи без аннотаций:

460486  461778  -
461863  463260  -
463562  464254  -
464413  465696  +
464563  465699  +
465782  466840  -
466914  467414  +

[464563 - 465699] — Вот нашёлся перекрывающийся участок на (+) цепи. Длина перекрывания — 1134 нуклеотида. Последовательность ORF длиннее на 147 нуклеотидов. Перекрывающаяся CDS кодирует продукт site-specific DNA-methyltransferase, то есть гомолог m4C-МТазы. Теперь уверенно можем сказать, что среди предсказанных ORF мы нашли CDS (WP_010700950.1), кодирующую гомолог метилтрансферазы. В NCBI у WP_010700950.1 не указан EC в аннотации.

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

Интересно поискать среди CDS генома последовательность, которая кодирует то же, что и наша ORF. Будем искать по аннотациям CDS.

elink отбирает UID (уникальные идентификаторы) белков, связанных с данной нуклеотидной записью через систему Entrez
efilter Фильтрует список белков, оставляя только те, у которых указан EC‑код: 2.1.1.37 — m5C‑метилтрансфераза; 2.1.1.72 — m6A‑метилтрансфераза; 2.1.1.113 — m4C‑метилтрансфераза.
efetch извлекает последовательности белков и выводит их ACCESSION 

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

После поиска по всем трём метилтрансферазам что-то нашлось только для m6A-MTазы.
Нашлось две записи, две ДНК-аденинметилазы: WP_022708611.1 (EC_number="2.1.1.72") и WP_010701371.1 (EC_number="2.1.1.72"). Ни та, ни другая не являются белком, который мы нашли в предыдущем пункте, но выполняют схожую функцию.