В прошлом семестре я работал с референсным протеомом бактерии Streptococcus gordonii str. Challis substr. CH1. Этому протеому UP000001131 соответствует код доступа GenBank GCA_000017005.1. Далее по этому коду был проведен поиск в базе данных NCBI Genomes. Была найдена соответствующая геномная сборка. Идентификатор RefSeq — GCF_000017005.1.
TaxID исследуемой бактерии можно узнать по полю OX в файле протеома. В нашем случае TaxID — 467705.
Чтобы скачать последовательности генома и таблицы локальных особенностей сборки GCF_000017005.1, была составлена следующая команда:
datasets download genome accession GCF_000017005.1 --include gff3,genome
Пояснения:
--include gff3,genome — получить файл с последовательностью генома и таблицу особенностей (feature table) в формате gff3
Архив был разархивирован с помощью следующей команды:
unzip ncbi_dataset.zip
Прежде чем проводить поиск открытых рамок считывания (ORF), нужно узнать, какой вариант генетического кода используется исследуемой бактерией. Для этого я скачал запись Streptococcus gordonii str. Challis substr. CH1 из базы данных NCBI Taxonomy. Для этого была использована следующая команда:
efetch -db 'taxonomy' -id '467705' -format 'xml'
Пояснения:
-db — база данных, в которых проводится поиск
-id — TaxID бактерии, который я узнал из файла генома по полю OX в файле протеома
-format — формат выходного файла (xml)
Информация о используемой бактерией генетической таблице находится в поле GeneticCode. Она использует 11 таблицу, что довольно предсказуемо, так как это самый распространенный вариант генетического кода для бактерий, архей и пластид.
Команда для получения последовательностей открытых рамок считывания:
getorf -sequence ncbi_dataset/data/GCF_000017005.1/GCF_000017005.1_ASM1700v1_genomic.fna -outseq ORFs.fasta -table 11 -minsize 150 -find 0
Пояснения:
-sequence — название файла генома бактерии (можно -sequence не указывать, я указал его для удобства восприятия)
-outseq — название выходного файла, в котором будут записаны все ORFs (тут также опцию можно не указывать)
-table — номер таблицы генетического кода, используемой бактерией (в моем случае 11)
-minsize 150 — выводить ORFs с не менее 50 аминокислотными остатками (150 нуклеотидов = 50 АК)
-find 0 — выводить ORFs между стоп-кодонами (в help написано: 0 Translation of regions between STOP codons); можно этот параметр не указывать, так как 0 стоит по умолчанию
Для того, чтобы проверить, что во всех ORFs находятся не менее 50 аминокислотных остатков, я составил следующий конвейер:
infoseq ORFs.fasta -only -length | sort -n | head
Пояснения:
-only -length — нужно вывести только длину последовательностей
sort -n сортирует длины по убыванию, head выводит только первые несколько строчек (больше и не надо). Минимальная длина оказалось равной 50, значит конвейер работает правильно.
Команда для создания базы белков для blastp:
makeblastdb -in ORFs.fasta -dbtype prot -out proteome
Пояснения
-in ORFs.fasta — название файла со всеми ORF длиной не менее 50 АК
-dbtype prot — указвыаем, что база данных именно белковая
-out proteome — имя базы данных
Конвейер для записывания в файл кодов доступа swissprot метилтрансфераз:
echo -e 'sw:P0AED9\nsw:P0AEE8\nsw:P23941' > methyltransferases.txt
Команда для выдачи последовательности этих метилтрансфераз, которые были записаны в файл methyltransferases.txt:
seqret @methyltransferases.txt query.fasta
Можно это сделать и при помощи одного конвейера:
echo -e 'sw:P0AED9\nsw:P0AEE8\nsw:P23941' | seqret -filter @stdin query.fasta
С помощью ранее составленной базы данных proteome и последоавтельностей methyltr_sequences.fasta была найдена ORF, гомологичная метилтрансферазе.
blastp -query query.fasta -db proteome -out methyltr_blastp.out -outfmt 7
Пояснения:
-query query.fasta — имя файла с последовательностями метилтрансфераз
-db proteome — ранее созданная мной база белков из ORF
-out methyltr_blastp.out — имя выходного файла
-outfmt 7 — табличная выдача в формате 7 (Tabular with comment lines)
Лучшая по весу (142) находка имеет идентификатор NC_009785.1_1552. Она гомологична m5C-МТазе E.coli, P0AED9.
Чтобы найти координаты NC_009785.1_1552 в геноме бактерии, я составил конвейер:
cat ORFs.fasta | grep 'NC_009785.1_1552'
Координаты: [689601-690914].
Далее я сделал поиск близких последовательностей по координатам. Для этого мне нужно было составить таблицу, которая содержит CDS той же последовательности, что и находка. Был использован конвейер:
cat ncbi_dataset/data/GCF_000017005.1/genomic.gff | grep 'NC_009785.1' | grep 'CDS' | cut -f4,5,7,9 > CDS_NC_009785_1.tsv
Затем был создан файл neighbors.tsv, в котором находятся близкие по координатам последовательности:
echo -e '689601\t690914\t-\tFOUND-ORF' | cat - CDS_NC_009785_1.tsv | sort -n | grep -C 3 'FOUND-ORF' > neighbors.tsv
Ссылка на файл neighbors.tsv
Из полученной таблицы видно, что была найдена последовательность, координаты которой практически полностью совпадают с находкой (689697-690917). Она принадлежит записи WP_006155197. В описании данной записи указано product=DNA cytosine methyltransferase, поэтому, скорее всего, эта последовательность действительно является гомологом.
Далее нужно было выяснить, смог ли я найти CDS, соответствующий моей находке, с помощью поиска по аннотации кодирующих участков в геноме. Для этого был составлен следующий конвейер:
elink -db nuccore -target protein -id 'NC_009785.1' | efilter -query 'EC[ECNO]' | esummary
Пояснения:
elink -db nuccore -target protein -id 'NC_009785.1' — ищет в базе данных Nucleotides запись с идентификатором NC_009785.1 и находит все соответствующе белковые последовательности в базе данных Protein.
efilter -query 'EC'[ECNO] — поиск соответствующих метилтрансфераз
esummary — 'синоним' efetch c опцией -format 'docsum', для более удобного просмотра полученных данных
В поле EC я указывал 2.1.1.37 (у m5C), 2.1.1.72 (у m6A) и 2.1.1.113 (у m4C). По m5C ничего не было найдено. По m6A была найдена 1 запись (WP_012000059.1), по m4C тоже 1 запись (WP_012000118.1). Среди этих записей моей находки не оказалось.