В этом практикуме я буду использовать протеом 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_23960, ее вес 628, координаты [2193705-2192257].
Далее определим с какими участками (CDS) пересекается лучшая находка, для этого создаем вспомогательный файл CDS.tsv, содержащий информацию о CDS, находящихся в той же геномной последовательности, что и лучшая находка:
grep 'NC_013961.1' genomic.gff | grep CDS | cut -f4,5,7,9 > CDS.tsv
Ну и находим участки, которые по координатам наиболее близки с моей находкой:
echo -e '2193705\t2192257\t-\tMY-ORF' | cat - CDS.tsv | sort -n | grep -C 3 'MY-ORF' > neighbors.tsv
Моей находке соответствует RefSeq:WP_004158132.1 DNA cytosine methyltransferase(m5C) с координатами [2192254-2193663].
Команды для поиска CDS в аннотации генома Erwinia amylovora. Для этого используем ЕС-код 2.1.1.72 (m6A), 2.1.1.37 (m5C), 2.1.1.113 (m4C).
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
elink -target protein -db nuccore -id NC_013961.1 | efilter -query '2.1.1.37' | efetch -format 'fasta'
Выдача: >WP_004154781.1 DNA cytosine methyltransferase [Erwinia amylovora]
MKFEGVNKPTVVSMFSGCGGMDLGFVQAGYDVVWANDFDADACLTYKRNIGDIIHGDVTTLDVPDVKNLD
VLTAGFPCQPFSNAGSRKGINDPRGQLYKQTFKFIKKLSPKVVVFENVRGLLSTKVENGKLIDEIVDTLT
NSLGYHVMYRLVNFSHFGVPQNRIRVILIAIKDLEYLPYIFPEVEMNKDLTLAKTLAGITDELPNQNELM
KLNPQALHYGSMIPAGGSWKNIPYELLPDRWKKIRDDMPRYHYPNFFRRYDLHEIQGTITAAFKPENAGV
WHPVEDRIYSVREIARFQTFPDDFVFEGRSIKSKYQQIGNAVPPLIGRKIAEQIKNYIHKVELSKIQRTM
VESTHLNVNQPIHVQKVSLALF
Этот белок имеет ту же функцию, но ID разные, значит не вышло найти CDS соответсвующий моей находке..
elink -target protein -db nuccore -id NC_013961.1 | efilter -query '2.1.1.113' | efetch -format 'fasta'
Выдачи нет.