практикум №9

EMBOSS, EDirect, Datasets CLI, blast+

Поиск ДНК-метилтрансферазы в геноме бактерии по последовательности и по аннотации


В данном практикуме использован геном бактерии Photobacterium ganghwense, протеом которой (UP000035909) изучался в прошлом семестре.

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

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

zcat ../../term2/pr8/UP000035909.swiss.gz | grep '^OX' | awk '{print $2}' | uniq -c

Выдача подтверждает, что все белки данного протеома аннотированы для одно и того же организма:

   4693 NCBI_TaxID=320778

Для получения списка AC геномных сборок, соотвествующих данному таксону, была использована команда:

datasets summary genome taxon 320778 --as-json-lines | dataformat tsv genome --fields 'accession,organism-name,assminfo-name,organism-infraspecific-strain'
 Assembly Accession      Organism Name   Assembly Name   Organism Infraspecific Names Strain
 GCF_017329545.1 Photobacterium ganghwense       ASM1732954v1    C2.2
 GCF_001029455.1 Photobacterium ganghwense       ASM102945v1     DSM 22954
 GCF_003025595.1 Photobacterium ganghwense       ASM302559v1     JCM 12487
 GCF_019084105.1 Photobacterium ganghwense       ASM1908410v1    K3B
 GCA_001029455.1 Photobacterium ganghwense       ASM102945v1     DSM 22954
 GCA_003025595.1 Photobacterium ganghwense       ASM302559v1     JCM 12487
 GCA_017329545.1 Photobacterium ganghwense       ASM1732954v1    C2.2
 GCA_019084105.1 Photobacterium ganghwense       ASM1908410v1    K3B

Выдача содержала 4 индентификатора геномной сборки для исследуемой бактерии. Каждая сборка имеет версию в GenBank и RefSeq. Поиск по полю Assembly Name показал, что это разные сборки, а не версии одной сборки. После добавления поля Organism Infraspecific Names Strain стало понятно, что это геномные сборки для разных штаммов. Далее будет использоваться сборка GCF_001029455.1, поскольку она принадлежит тому же штамму, что и протеом (DSM 22954), и относится к RefSeq, то есть лучше аннотирована.

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

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

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

Полученный архив был распакован командой unzip.

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

Исследуемая бактерия "использует" таблицу генетического кода №11, что было выяснено командой:

efetch -db 'taxonomy' -id '320778' -format 'xml' | grep -A '1' 'GeneticCode'

Поиск и трансляция открытых рамок считывания между стоп-кодонами были осуществлены командой:

getorf GCF_001029455.1_ASM102945v1_genomic.fna -outseq translated_orfs.fasta -table 11 -minsize 150 -find 0

Пояснения к команде:

  • GCF_001029455.1_ASM102945v1_genomic.fna - фаста-файл с нуклеотидной п-тью генома;
  • translated_orfs.fasta - фаста-файл с аминокислотной п-тью транслированных ORF;
  • -table 11 - указание таблицы ген. кода;
  • -minsize 150 - указание минимальной длины ORF (для получения амк. п-тей, не короче 50 амк);
  • -find 0 - существует два возможных определения открытой рамки считывания: это может быть либо область, не содержащая стоп-кодонов, либо область, начинающаяся со стартового кодона и заканчивающаяся стоп-кодоном (из getorf -help). Соответственно, значение 0 - это транслирование п-ти между стоп-кодонами.

Для проверки того, что транслированные ORF действительно не превышают 50 амк по длине (так и получилось), была запущена программа:

infoseq translated_orfs.fasta -only -length | sort -n | less

Белковая база для алгоритма blastp на основе п-тей трансляций была создана командой:

makeblastdb -in translated_orfs.fasta -dbtype 'prot'

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

Поиск гомологов ДНК-метилтрансферазы исследуемой бактерии посредством белок-белкового бласта проводился среди таких бактерий, как E.coli и Bacillus amyloliquefaciens. Последовательности белков одним фаста-файлом были получены командой seqret из базы данных Swiss-Prot по следующим кодам доступа и командой:

для E.coli

для Bacillus amyloliquefaciens

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

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

Поиск гомологов метилтрансфераз по базе транслированных ORF генома исследуемой бактерии локальным алгоритмом blastp, где в качестве запроса дан список указанных выше белковых п-тей, был выполнен командой:

blastp -task blastp -query query_MTases.fasta -db translated_orfs.fasta -out blast_table.out -outfmt 7

С результатом выдачи можно ознакомиться по ссылке. Лучшая по весу и e-value находка имеет идентификатор NZ_LDOU01000024.1_2100.

Таблица 1. Характеристики лучшего выравнивания.
Название рамки Координаты в геноме Вес находки Гомолог МТазы...
NZ_LDOU01000024.1_2100 [28678 - 27854] (REVERSE SENSE) 357 m6A-МТаза E.coli

Координаты в геноме были получены командой grep по названию рамки в файле с транслированными рамками.

Чтобы определить, какие CDS из таблицы локальных особенностей генома располагаются рядом с лучшей находкой, сперва необходимо было получить строки из таблицы, которые соответствуют CDS в той же геномной последовательности, и столбцы с координатами (4 и 5), цепью (7) и дополнительной информацией (9). Данная информация была записана в файл CDS.tsv

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

Далее в данный файл была добавлена строка с координатами находки (с учетом обратных координат). С помощью сортировки получилось определить, какие 3 рамки находятся в окрестностях (с обеих сторон).

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

С результатом выдачи можно ознакомиться по ссылке. По результатам выдачи видно, что анализируемая находка почти полностью соответствует CDS с координатами 27851-28678, аннотированной как аденин-специфичная ДНК-метилтрансфераза, то есть найден гомолог m6A-МТазы. С остальными CDS находка не пересекается.

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

Поиск соответствующих находке CDS по аннотациям в геноме проводился по EC-коду подходящей метилтрансферазы (2.1.1.72 для m6A-МТазы). Поиск был выполнен командами edirect:

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

В выдаче был только один белок, совпадающий с таковым в предыдущем пункте (то есть имеющий тот же CDS id и такую же аннтоцию каталитической активности). Ниже приведена выдача без п-ти. Таким образом, гомолог метилтрансферазы у исследуемой бактерии удалось найти как по транслированному протеому, так и по геному. 

>WP_047887153.1 Dam family site-specific DNA-(adenine-N6)-methyltransferase [Photobacterium ganghwense]

Поиск по ЕС других метилтрансфераз из п.4 не дал результатов.