Задача практикума - найти одну из ДНК-метилтрансфераз в геноме бактерий или архей, основываясь на последовательности и аннотации и используя инструменты EMBOSS, EDirect, NCBI Datasets CLI и BLAST+
В прошлом семестре я работала с Pseudoalteromonas piscicida, тогда Proteome ID был UP000016521, на данный момент он стал избыточным(Redundant proteome). Поэтому я решила взять UP000305423 (я пыталась работать с UP000016521, но ничего не получалось, файл пустой)
AC сборки нам нужен для последующей загрузки последлвательности и таблицы локальных особенностей
На первом этапе для получения значения поля TaxID испоьзовалась следующая команда:
zcat UP000305423.swiss.gz | grep "^OX"| cut -c 17-22|uniq
Для получения списка сборок я испоьзовала esearch+efetch:
esearch -db assembly -query "txid43662[Organism]" | efetch -format docsum | xtract -pattern DocumentSummary -element AssemblyAccession,AssemblyName,AssemblyStatus > eta1.txt
Я выбрала сборку GCF_019797925.1 (ASM1979792v1) она одна референсная
Последовательность была скачана с помощью команды:
datasets download genome accession GCF_019797925.1 --include 'gff3','genome'
Далее датасет был распакован с помощью unzip
Чтобы узнать какой генетический код использует организм я использовала NCBI Тaxonomy
efetch -db taxonomy -id 43662 -format xml > ee.xml
Pseudoalteromonas piscicida использует таблицу 11
Открытые рамки считывание искалиссь с помощью команды:
getorf -sequence ncbi_dataset/data/GCF_019797925.1/GCF_019797925.1_ASM1979792v1_genomic.fna -outseq orfs.fasta -minsize 150 -find 0 -table 11
Чтобы проверить правильно ли указана длина рамок считывания исползовалась программа infoseq:
infoseq -only -length orfs.fasta| sort -n| head -n 2
Далее мы создаем белковую базу данных:
makeblastdb -in orfs.fasta -dbtype prot -out ORFs
Для скачивания последовательности ДНК-метилтрансферазы из Swiss-Prot использовали следующую команду:
echo -e "sw:P0AED9\nsw:P0AEE8\nsw:P23941" | seqret -sequence @stdin -outseq query_MTases.fasta
С помощью ниже преведенной команды выявилась находка, имеющая вес 328 (с названием:NZ_CP081860.1_1062, координатами [424868 - 425749], гомологом m6A)
blastp -query query_MTases.fasta -db ORFs -outfmt 7 -out results.txt
По координатам находок можно определить какие CDS из таблицы локальных особенностей будут рядом. Отбирали нужные столбцы с помощью команды:
grep 'NZ_CP081860.1' genomic.gff| grep 'CDS'| cut -f 4,5,7,9 > CDS.tsv
Отбирали соседние CDS ткаим образом:
echo -e '424868\t425749\t-\tMY-ORF' | cat - CDS.tsv | sort -n | grep -C 3 'MY-ORF' > neighbors.tsv
Нашлось пересечение c ID=cds-WP_010368897.1 (координаты CDS 424919-425752, координаты пересечения:424919-425749)
Для того, чтоб узнать, можно ли найти CDS, соответсвующей найденной находке использовался конвейер:
elink -target protein -db nuccore -id "NZ_CP081860.1" | efilter -query "(2.1.1.37 OR 2.1.1.72 OR 2.1.1.113)" | efetch -format docsum > hehe.txt
Был найден 1 белок WP_010368897.1,это как раз тот кодирующий участок, для которого найдены пересечения