Для анализа был выбран протеом археи Methanosarcina vacuolata Z-761(UP000033096). Её TaxID был получен следующей командой:
grep 'TaxID' UP000033096.swiss | cut -b1-24 | uniq -c
Выдача была следующей:
3478 OX NCBI_TaxID=1434123
Поиск AC генома осуществлялся с помощью следующей команды;
datasets summary genome taxon 1434123 --as-json-lines | dataformat tsv genome | vd
Всего была найдена одна сборка генома в двух вариантах: GenBank и RefSeq. Была выбрана сборка RefSeq(GCF_000969905.1).
Скачивание геномной сборки и распаковка скачанного архива осуществлялась с помощью следующих команд:
datasets download genome accession GCF_000969905.1 --include 'genome' --include 'gff3'
unzip ncbi_dataset.zip
Далее требуется выяснить какой вариант генетического кода использует выбранная архея. Для этого использовалась команда:
esearch -db taxonomy -query '1434123[UID]' | efetch -format native -mode xml | less
Оказалось, что архея использует 11-й вариант генетического кода(Bacterial, Archaeal and Plant Plastid). Получение открытых рамок считывания от стоп-кодона до стоп-кодона осуществлялось с помощью команды:
getorf GCF_000969905.1_ASM96990v1_genomic.fna -outseq proteins.fasta -find 0 -minsize 150 -table 11
Белкова база данных по открытым рамкам считывания была создана при помощи команды:
makeblastdb -in proteins.fasta -dbtype 'prot' -out ORFs
Далее было проверено, действительно ли среди полученных рамок считывания нет тех, длина которых меньше 50 аминокислот:
infoseq proteins.fasta -filter | cut -f4 -dP|cut -b5-15 |sort -n|less
Сначала с помощью текстового редактора nano был создан файл(mt.txt), содержащий USA всех трех последовательностей, записанных по строкам(sw:P0AED9, sw:P0AEE8, sw:P23941). Далее эти последовательности были скачаны с помощью команды:
seqret @mt.txt -filter > query_MTases.fasta
Поиск гомологов метилтрансфераз осуществлялся при помощи blastp в созданной базе ORFs следующей командой:
blastp -query query_MTases.fasta -out out.txt -db ORFs -outfmt 7
Результаты поиска доступны по ссылке.
У лучшей находки гомолога M6A-метилтрансферазы был идентификатор NZ_CP009520.1_20066, вес - 46.2 и координаты в геноме - [3153279-3152356](REVERSE SENSE). Теперь найдем CDS в геноме с похожими координатами на нашу находку. Сперва, выделим из файла локальных особенностей информацию о координатах и аннотацию всех CDS командой:
cat genomic.gff | grep 'CDS' | cut -f4,5,7,9 > CDS.tsv
Поиск CDS соседних или перекрывающихся с моей находкой осуществялся с помощью следующей команды:
echo -e '3152356\t3153279\t-\tMY-ORF' | cat - CDS.tsv | sort -n | grep -C 3 'MY-ORF' > neighbors.tsv
Файл выдачи можно посмотреть по ссылке.
Результаты: найденная открытая рамка считывания(3152356 3153279 (-) ) пересекается с CDS(3152353 3153258 (-) ) гена, кодирующего ДНК-аденинметилтранферазу.
Поиск гомологов метилтрансфераз осуществлялся с помощью следующих трех команд:
esearch -db nuccore -query 'GCF_000969905.1' | elink -target protein |efilter -query '2.1.1.37' | efetch -format fasta | less
esearch -db nuccore -query 'GCF_000969905.1' | elink -target protein |efilter -query '2.1.1.72' | efetch -format fasta | less
esearch -db nuccore -query 'GCF_000969905.1' | elink -target protein |efilter -query '2.1.1.113' | efetch -format fasta | less
Результаты: ни одного белка по данным трем запросам найдено не было.