В этом практикуме я использовал геном археи Thermococcus kodakarensis KOD1 (NCBI Taxonomy ID 69014). Ее протеом уже был скачан в формате swiss; чтобы найти в скачанном файле TaxID организма, а также убедиться, что его значения совпадают у всех белков, я использовал следующий скрипт:
zcat UP000000536.swiss.gz | grep 'TaxID' | cut -c 1-21 | uniq
Выдача скрипта:
OX NCBI_TaxID=69014
С помощью программы datasets я получил список сборок, соответствующих T. kodakarensis KOD1 в RefSeq, для этого был использован следующий скрипт:
datasets summary genome taxon 69014 --assembly-source refseq --as-json-lines | dataformat tsv genome --fields accession
Скрипт выдал единственную сборку GCF_000009965.1. Команда для скачивания последовательности генома и таблицы локадьных особенностей:
datasets download genome accession GCF_000009965.1 --include genome,gff3
Полученный архив был распакован с помощью следующей команды:
unzip ncbi_dataset.zip
Для дальнейшей работы с геномом T. kodakarensis нужно было проверить, какой он использует генетический код. Это было сделано с помощью команды efetch:
efetch -db taxonomy -id 69014 -format xml > taxon.xml
Из созданного xml-документа видно, что таблица кода 11. С помощью программы getorf были получены трансляции открытых рамок считывания между стоп-кодонами длиной не менее 50 аминокислот:
getorf ncbi_dataset/data/GCF_000009965.1/GCF_000009965.1_ASM996v1_genomic.fna seq.fasta -find 0 -minsize 150 -table 11
По промежуточному файлу seq.fasta была создана база данных ORFs:
makeblastdb -dbtype prot -in seq.fasta -out ORFs
Чтобы убедиться, что все трансляции имеют длину не менее 50 аминокислот я воспользовался программой infoseq:
infoseq seq.fasta |cut -d ' ' -f35 | sort -n | uniq
С помощью команды seqret были скачаны последовательности трех метилтрансфераз из Swiss-Prot, для этого понадобилось создать также промежуточный файл listfile.txt:
echo -e 'sw:P0AED9\nsw:P0AEE8\nsw:P23941' > listfile.txt
seqret @listfile.txt query_MTases.fasta
Я провел выравнивание последовательностей из базы данных ORFs с последовательностями из query_MTases.fasta с помощью локального blastp:
blastp -db ORFs -query query_MTases.fasta -out blast.txt -outfmt 7
Наибольший вес (39.3) имеет выравнивание рамки NC_006624.1_7739 (координаты - 1835582 - 1836472) с m4C-метилтрансферазой Bacillus amyloliquefaciens.
Чтобы найти CDS, пересекающиеся с выравненной рамкой считывания, я создал таблицу, содержащую строки из таблицы локальных особенностей T. kodakarensis, соответствующие геномной последовательности NC_006624.1, и столбцы с координатами, цепью и дополнительной информацией:
grep NC_006624.1 ncbi_dataset/data/GCF_000009965.1/genomic.gff | grep CDS | cut -f4,5,7,9 > CDS.tsv
Затем я вставил в таблицу свою ORF и отобрал строки с самыми близкими CDS:
echo -e '1835582\t1836472\t-\tMY-ORF' | cat - CDS.tsv | sort -n | grep -C 3 'MY-ORF' > neighbors.tsv
Из получившейся таблицы видно, что NC_006624.1_7739 почти полностью перекрывается с геном белка, содержащего метилтрансферазный домен (cds-WP_011250995.1; координаты - 1835603 - 1836475).
Я также проверил, возможно ли найти CDS, соответствующие NC_006624.1_7739, с помощью поиска по аннотации кодирующих участков в геноме. Для этого был использован конвейер edirect:
elink -db nuccore -target protein -id NC_006624.1 | efilter -query '2.1.1.113' | efetch -format 'fasta'
2.1.1.113 - m4C-метилтрансфераза B. amyloliquefaciens, кроме нее, я искал по аннотациям m5C-метилтрансферазы E. coli (2.1.1.37) и m6A-метилтрансферазы E. coli (2.1.1.72). Ни один из вариантов поиска не дал результатов.