В этом практикуме я буду использовать протеом UP000001841 бактерии Erwinia amylovora (strain CFBP1430), так как с ним я уже работала в практикуме 8 прошлого семестра и эти файлы у меня уже есть.
Сначала я получила TaxID организма, которому принадлежит геном, распаковав файл и применив команду:
grep '^OX' UP000001841.swiss | cut -d ' ' -f4 | tr -d 'NCBI_TaxID=;' | sort -u
Выдача: 665029
Получаем список AC сборок необходимого генома с помощью команды:
esearch -db assembly -query '665029' | efetch -format 'docsum' | xtract -pattern DocumentSummary -element AssemblyAccession,TaxId,Title
Выдача: GCF_000091565.1
Геномная сборка всего одна, ее и используем
Имея AC геномной сборки, скачаем последовательность и feature table с помощью команды:
datasets download genome accession GCF_000091565.1 --include genome,gff3 --filename feature_table.zip
И распакуем файл с помощью unzip
Прежде, чем искать открытые рамки считывания, нужно определить, какой вариант генетического кода использует исследуемый организм. Сделаем это с помощью команды:
efetch -db taxonomy -id '665029' -format 'xml' | grep -H '
Выдача: GCId 11 GCId
Теперь мы знаем, что рассматриваемая бактерия использует таблицу 11.
Ищем открытые рамки считывания и получаем их трансляции с помощью программы getorf из пакета EMBOSS.
getorf -sequence 'GCF_000091565.1_ASM9156v1_genomic.fna' orfs.fasta -minsize 150 -table 11 -filter
Проверим, что их трансляции не короче 50 аминокислот.
infoseq -sequence orfs.fasta -only -length | sort -n -u
Создаем по ним белковую базу для blastp ORFs.
makeblastdb -in orfs.fasta -dbtype 'prot'
Скачать последовательности белков P0AED9 (Dcm, m5C-МТаза, E.coli), P0AEE8 (Dam, m6A-МТаза, E.coli), и P23941 (m4C-МТаза, Bacillus amyloliquefaciens) из Swiss-Prot можно с помощью команды:
echo -e 'sw:P0AED9\nsw:P0AEE8\nsw:P23941' | seqret -filter '@stdin' 'query_MTases.fasta'
С помощью локального blastp найдем белковые последовательности в геноме (ORFs), схожие с последовательностями ДНК-метилтрансфераз (query_MTases.fasta). Команда:
blastp -query query_MTases.fasta -db orfs.fasta -out orfsblastp.out -outfmt 7
Теперь найдем находку с наибольшим весом: NC_013961.1_18395, ее вес 407, координаты [3521845 - 3521012], гомолог m6A-МТазы.
Далее определим с какими участками (CDS) пересекается лучшая находка, для этого создаем вспомогательный файл CDS.tsv, содержащий информацию о CDS, находящихся в той же геномной последовательности, что и лучшая находка:
grep 'NC_013961.1' genomic.gff | grep CDS | cut -f4,5,7,9 > CDS.tsv
Ну и находим участки, которые по координатам наиболее близки с моей находкой:
echo -e '3521845\t3521012\t-\tMY-ORF' | cat - CDS.tsv | sort -n | grep -C 3 'MY-ORF' > neighbors.tsv
Из файла можно увидеть, что моя находка соответствует CDS site-specific DNA-methyltransferase (adenine-specific) с координатами [3521009-3521830], то есть m6A-МТазе.
В завершении, найдем CDS, соответствующие моей находке, с помощью поиска по аннотации кодирующих участков в геноме. Для этого используем ЕС-код 2.1.1.72 (m6A), так как в прошлом задании мы выяснили, что именно он соответствует находке, с другими EC никаких результатов не будет.
elink -target protein -db nuccore -id NC_013961.1 | efilter -query '2.1.1.72' | efetch -format 'fasta'
Выдача: >WP_004160689.1 adenine-specific DNA-methyltransferase [Erwinia amylovora]
MKKNRAFLKWAGGKYPLLDDIRSHLPEGDCLIEPFVGAGSVFLNTQYSRYILADINSDLINLYNIVKERT
DEFVLDARLLFTPESNDADVYYAWRQEFNQCGDAYRRALLFLYLNRHCYNGLCRYNLSGEFNVPFGRYRR
PYFPQEELYWFAERAQYATFVCESYDVTLGKAHQGSVVYCDPPYAPLSATANFTAYHTNSFSMREQQHLA
ELATRLAQESSITVLISNHDTPLTRQWYQGATKLQEIKVRRSISRSGNSRSKVNELLALYNGR
ВСЕ!