Практикум №9. EMBOSS, EDirect, NCBI Datasets CLI и blast+

Введение

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

Непосредственной задачей был поиск ДНК-метилтрансферазы в геноме бактерии Streptomyces collinus по её последовательности и аннотации.

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

Перед началом работы используя идентификатор протеома: UP000015423 на странице базы данных UniProt Proteomes, а также базы NCBI Assembly (Datasets Genome) были получены следующие данные:

Код доступа геномной сборки из GenBank (GCA_000444875.1)

Страница сборки в Genome

Код доступа из RefSeq (GCF_000444875.1)

Последовательность генома и таблица локальных особенностей

По AC геномной сборки RefSeq был скачан геном и таблица локальных особенностей в формате GFF3.

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

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

Перед дальнейшей работой нужно было определить вариант используемого генетического кода исследуемой бактерией. Для этого с помощью efetch была скачана запись про таксон в базе NCBI Тaxonomy в формате xml:

efetch -db taxonomy -id 1214242 -format xml > output.xml

id был взят со страницы таксона на сайте NCBI. В абзаце GeneticCode было указано значение 11, что соответствует стандартной таблице №11.

Далее используя программу getorf из пакета EMBOSS были найдены открытые рамки считывания и сразу получены их трансляции:

getorf -sequence GCF_000444875.1_ASM44487v1_genomic.fna -table 11 -minsize 150 -outseq peptides

Из полученного списка белков была создана база данных proteome:

makeblastdb -in peptides -dbtype prot -out proteome

Для проверки порога длины рамок считывания из EMBOSS была запущена программа infoseq:

infoseq -sequence peptides -only -length -noheading | sort -nr | tail -n 1 > output.txt

Самая короткая последовательность составляет 50 а.о., что не противоречит условию.

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

При поиске сходства последовательностей метилтрансфераз в качестве сравниваемого фермента было высказано предположение, что для m4C-МТаза бактерии Bacillus amyloliquefaciens с большей вероятностью обнаружится сходство в последовательностях ферментов. Данный выбор был обоснован большим родством организмов, чем с E.coli.

Для загрузки одним конвеером использовался следующий запрос:

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

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

При помощи программы blastp произведён поиск в полученной базе данных схожих с метилтрансферазой последовательности среди белков иследуемой бактерии:

blastp -query query.fasta -db proteome -out results.txt -outfmt 7

Полученный файл results.txt

Если оценивать находки по проценту идентичности, тогда лучшей из них является ДНК-цитозин метилтрансфераза E.coli. Однако стоит принять во внимание значения evalue, которые для всех последовательностей являются слишком большими.

Из низкой статистической значимости сходства двух сравниваемых последовательностей следует высокая вероятность случайного появления наблюдаемого уровня сходства. Из этого можно сделать вывод, что последовательности базы данных и схожих с метилтрансферазой последовательности среди белков иследуемой бактерии нельзя считать гомологичными.

В ходе дальнейшей работы были установлены координаты "лучшей" находки:

grep 'NC_021985.1_49758' peptides

А также получена таблица CDS:

grep 'CDS' genomic.gff | cut -f 4,5,7,9 > CDS.tsv

По координатам (6085197 - 6082249) был произведен поиск расположенных рядом CDS

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

Полученный файл neighbors.tsv

Из выдачи следует, что наиболее перекрывающимся геном является белок-повтор тетратрикопептида системы FxSxx-COOH с координатами (6082414-6085242).

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

Так же в работе было проверена возможность нахождения "лучшей" находки по аннотированным генам исследуемой бактерии, для чего был сделан запрос entrez:

elink -db nuccore -id 'NC_021985.1' -target 'protein' | efilter -query '2.1.1.37[EC/RN Number]' | efetch -format 'ipg' | less

По пустой выдаче можно судить, что найти данный участок по аннотации не представляется возможным.