Практикум 9. EMBOSS, Entrez Direct, NCBI Datasets

Протеом бактерии Thiomicrorhabdus sediminis (G1) (ID штамма: 2580412) сохранился с прошлого семестра (Практикум 8). В дальнейшем я буду работать с ним (и списывать у самого себя из прошлого). Идентификатор протеома: UP000304864

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

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

В UniProt Proteomes я выполнил поиск протеома по его идентификатору UP000304864 и на странице записи о нем узнал Taxon ID организма: 2580412, а также ID сборки на RefSeq: GCA_005885815.1. По этому коду был выполнен поиск в NCBI Assembly, где была найдена соответствующая сборка в RefSeq с идентификатором GCF_005885815.1. В дальнейшем будет использоваться именно эта сборка из RefSeq.

Идентификаторы и ссылки в явном виде

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

Использованные команды:
datasets download genome accession GCF_005885815.1 --include gff3,genome
unzip ncbi_dataset.zip
Для удобства работа велась в папке public_html/term3/pr9/dataset/.

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

С помощью команд
efetch -db 'taxonomy' -id '2580412' -format 'xml' > entry.xml
cat entry.xml | grep 'GCId'
было получено, что Thiomicrorhabdus sediminis использует 11 таблицу генетического кода.

Поиск и трансляция открытых рамок считывания (Open Reading Frame, ORF) были выполнены с помощью команды
getorf GCF_005885815.1_ASM588581v1_genomic.fna -outseq translated_orfs.fasta -table 11 -minsize 150 -find 0
Объяснение команды:

С помощью команды
infoseq translated_orfs.fasta -only -length | sort -n | less
было подтверждено, что транслированные ORF имеют длину не менее 50 аминокислот.

Белковая база для будущего blastp была создана с помощью команды
makeblastdb -in translated_orfs.fasta -out proteome -dbtype 'prot'

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

Я создал файл query.fasta, содержащий последовательности трех ДНК-метилтрансфераз: P0AED9 (Dcm, m5C-МТаза, E.coli), P0AEE8 (Dam, m6A-МТаза, E.coli), и P23941 (m4C-МТаза, Bacillus amyloliquefaciens). Использованная команда:
echo "sw:P0AED9 sw:P0AEE8 sw:P23941" | tr ' ' '\n' | seqret -filter '@stdin' -outseq query.fasta

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

С помощью метода blastp по локальной базе данных proteome был проведен поиск последовательностей, обладающих наибольшим сходством.
blastp -query query.fasta -outfmt 7 -db proteome -out blastp_out.txt
Файл с выдачей blastp: blastp_out.txt

Рамка NZ_CP040602.1_1456 имела наибольший вес – 75.9 (evalue: 9.43e-16). Обнаруженный гомолог: Dcm, m5C-МТаза, E.coli. Координаты: 460215 - 461321. Команда, использованная для поиска координат:
grep 'NZ_CP040602.1_1456 ' translated_orfs.fasta
Пробел после запроа в grep использовался для отсеивания результатов вида *_14560, *_14561 и т.п.

Из таблицы локальных особенностей были отобраны строки, соотвествующие CDS, и столбцы с координатами (4,5), цепью (7) и дополнительной информацией (9)
grep 'CDS' genomic.gff | cut -f1,3,4,5,7,9 > CDS.tsv
Комментарий: проверка длины трансляций была убрана из текущего пункта за ненадобностью.

Для поиска пересечений я использовал конвейер 

echo -e 'NZ_CP040602.1\tCDS\t460215\t461321\t+\tFOUND-ORF' | cat - CDS.tsv | sort -k3,3 -n | grep -C 3 'FOUND-ORF' > neighbors.tsv
В результате работы конвейера был создан файл neighbors.tsv. С FOUND-ORF пересекалась одна CDS с практически совпадающими координатами (460221–461324), у которой был указан продукт ДНК цитозин-метилтрансфераза (DNA cytosine methyltransferase). Небольшое отличие в координатах объясняется тем, что FOUND-ORF получена программой getorf с опцией поиска между стоп-кодонами (-find 0), в то время как аннотированная CDS предсказана с использованием более строгих правил: она должна начинаться со старт-кодона и заканчиваться стоп-кодоном. Разница в 6 нуклеотидов на старте (2 аминокислоты) укладывается в это объяснение: getorf включил в предсказание лидерную последовательность до старт-кодона. На 3'-конце отличие составляет 3 нуклеотида, что связано с тем, что getorf включает стоп-кодон в состав ORF (поскольку ищет фрагменты между стоп-кодонами), тогда как в аннотированной CDS стоп-кодон не входит в кодирующую последовательность. Таким образом, разница в координатах обусловлена особенностями алгоритмов предсказания, а не ошибкой аннотации.

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

Я выполнил поиск по базе данных nuccore с помощью слеующего конвейера: 

elink -db nuccore -id 'NZ_CP040602.1' -target protein | efilter -query "2.1.1.37 [ECNO]" | efetch -format acc

Пояснение:

Был найден белок WP_138563936.1, что означает совпадение результата с предыдущим.