Будем искать одну из ДНК-метилтрансфераз в геноме Thermococcus profundus, чья частичная характеристика приведена ниже
код доступа из RefSeq: GCF_002214585.1
код доступа сборки GenBank: GCA_002214585.1
ссылка на страницу сборки Genome.
Идентификатор протеома: UP000250179.
Скачивание данных. Эта команда загрузит геномную последовательность (в формате FASTA) и аннотации (в формате GFF3) для указанной сборки.
datasets download genome accession GCF_002214585.1 --include gff3,genome
--include gff3,genome — опция, указывающая, что нужно включить в архив файлы с аннотацией (GFF3) и геномной последовательностью.
в текущей директории появится файл с именем ncbi_dataset.zip. Его нужно распаковать.
unzip ncbi_dataset.zip
Первым делом определим, какой вариант генетического кода использует исследуемый организм. Команда для получения записи про таксон:
efetch -db taxonomy -id 49899 -format xml
-db — база данных, в которых проводится поиск
-id — TaxID бактерии cо страницы в Proteomes
-format — формат выходного файла (xml)
В поле GeneticCode найдем, что исследуемым организмом используется 11 вариант генетического кода.
Теперь используем программу getorf для поиска всех открытых рамок считывания в геномной последовательности.
getorf -sequence ncbi_dataset/data/GCF_002214585.1/GCF_002214585.1_ASM221458v1_genomic.fna -table 11 -minsize 150 -find 0 -outseq proteome_orf.faa
-minsize — минимальная длина ORF в нуклеотидах. 150 — соответствует 50 аминокислотам (150 ÷ 3 = 50). Отсекает случайные короткие ORF.
-find — режим поиска ORF:
0 —искать между стоп-кодонами
Создание белковой базы для BLAST
Файл proteome_orf.faa уже является fasta-файлом с белковыми последовательностями. Чтобы использовать его как базу для blastp, ее нужно отформатировать с помощью makeblastdb.
makeblastdb -in proteome_orf.faa -dbtype prot -out proteome
проверим, что в нашем файле proteome_orf.faa нет последовательностей короче 50 аминокислот.
infoseq proteome_orf.faa -only -length | sort -n | head
-only -length указывает, что нужно выводить только длины.
Для поиска гомологов были выбраны три ДНК-метилтрансферазы из базы данных Swiss-Prot.
Команда для скачивания всех трех последовательностей одним конвейером:
echo -e 'sw:P0AED9\nsw:P0AEE8\nsw:P23941' | seqret -filter @stdin -outseq query.fasta
echo -e - формирует список идентификаторов белков с переносами строк (\n)
sw: - префикс, указывающий на базу данных Swiss-Prot
seqret - программа для извлечения последовательностей из баз данных
-filter @stdin - чтение списка идентификаторов из стандартного ввода
-outseq query.fasta - сохранение результатов в файл
Найдём гомологи известных ДНК-метилтрансфераз в предсказанном протеоме Thermococcus profundus и проанализируем геномное окружение найденного гена.
blastp -query query.fasta -db proteome -out blast_results.txt -outfmt 7
Лучщая находка orf NZ_CP014862.1_153 с координатами 33493 - 34626, которые посмотрели командой
grep ">NZ_CP014862.1_153 " proteome_orf.faa
с весом 41,6 является гомологом m4C-метилтрансферазы (M.BamHI)
определим, какие CDS из таблицы локальных особенностей генома располагаются рядом.
отберем нужные строки и столбцы из таблицы локальных особенностей
grep '^NZ_CP014862.1' ncbi_dataset/data/GCF_002214585.1/*.gff | grep 'CDS' | cut -f 4,5,7,9 > CDS.tsv
Ищем близкие CDS
echo -e '33493\t34626\t+\tFOUND-ORF' | cat - CDS.tsv | sort -n | grep -C 3 'FOUND-ORF' > neighbors.tsv
ORF (33493-34626) пересекается с аннотированным геном WP_232473543.1 (33754-34629). В описании указана метилтрансферазная активность, поэтому можно предположить что последовательности гомологичны.
можно ли было найти нашу метилтрансферазу через поиск по EC-кодам в аннотированных генах? сравним с результатами BLAST-поиска.
elink -db nuccore -id NZ_CP014862.1 -target protein | efilter -query " 2.1.1.113" | efetch -format acc
elink -db nuccore -id NZ_CP014862.1 -target protein - получает все белковые записи, связанные с нашим геномом
efetch -format acc - получает accession numbers найденных белков
По результатам этих запросов не было найдено ни одной записи.