В данном практикуме я искала ДНК-метилтрансферазу в геноме Mycobacterium tuberculosis, упомянутой в практикуме 8 из 2 семестра (файл UP000001584.swiss).
Поиск осуществлялся с помощью с помощью EMBOSS, EDirect, NCBI Datasets CLI и blastp.
Для начала я получила ID такснона моей бактерии :
grep '^OX' UP000001584.swiss | tr ' ' '\t' | cut -f4 | tr -d ';' |sort | uniq
Затем с помощью этого ID нашла сборки генома
esearch -db 'taxonomy' -query ''83332' [UID]' | elink -target 'assembly' | efetch -format 'uid' - получение номеров сборок
datasets summary genome taxon 83332 --as-json-lines | dataformat tsv genome | cut -f1 | sort | uniq | less - а здесь получение АС сборок
(GCF_000195955.2 - нужная сборка, с которой я работаю в дальнейшем)
Зная АС сборки можно загрузить нужные нам таблицы :
datasets download genome accession GCF_000195955.2 --include gff3 --include genome
unzip ncbi_dataset.zip - распаковка таблицы, т.к она скачалась в формате архива.
Для начала, я определила какую таблицу генетического кода использует мой организм :
efetch -db taxonomy -id '83332' -format 'xml' | less - узнаем, что нам нужна 11 таблица
(GeneticCode 11)
Теперь можно получать открытые рамки считывания :
getorf 'GCF_000195955.2_ASM19595v2_genomic.fna' -table 11 -minsize 150 getorf.fasta - получения рамок считывания в getorf.fasta
Далее создадим базу данных blast...
makeblastdb -in getorf.fasta -out ORFs -dbtype prot
и проверим, что по длине они не меньше 50 нуклеотидов
infoseq -sequence getorf.fasta -only -length | sort -n| less (да, все рамки достаточно длинные)
Как нам известно, у всех ДНК-метилтрансферазы прокариот есть гомологичные каталитические домены. Для поиска ДНК-метилтрансферазы у нашего организма скачаем последовательности белков гомологов (чтобы затем осуществить поиск по генам) :
echo 'sw:P0AED9' 'sw:P0AEE8' 'sw:P23941'| tr ' ' '\n' > seq.txt
seqret @seq.txt query_MTases.fasta
Таким образом мы получили m5C-МТазу E.coli, m6A-МТазу E.coli и m4C-МТазу Bacillus amyloliquefaciens (поиск производился по их кодам доступа Swiss-Prot).
Применим blastp по ранее созданной базе данных ORFs и посмотрим на выдачу (нужна находка с наибольшим весом):
blastp -task blastp -query query_MTases.fasta -db ORFs -out blastp.out -outfmt 7
Возьмем CDS самой лучшей по весу (30.0) находки из файла blastp.out (это гомолог m5C-МТазы):
grep AL123456.3 genomic.gff |grep CDS| cut -f4,5,7,9 > CDS.tsv
В файле getorf.fasta находим координаты AL123456.3_33444 (33444 - номер рамки в геноме, AL123456.3 - AC находки), координаты [3141314 - 3142525]
Теперь ищем CDS, пересекающиеся с находкой, по их описанию можно понять, за что данный ген отвечает (в идеале, должна быть DNA-cytosine methyltransferase где-то в описании):
echo -e '3142525\t3141314\t-\tMY-ORF' | cat - CDS.tsv | sort -n | grep -C 3 'MY-ORF' > neighbors.tsv
Увы, ничего сильно полезного из выдачи файла neighbors.tsv я не получила. Находка практически полностью совпала по координатам с CCP45636.1
(координаты 3141311-314222), но, к сожалению, product этого CDS ничем интересным не отличается (это Probable Sn-glycerol-3-phosphate transport integral membrane protein ABC transporter UgpA - у него транспортная функция (перемещение субстрата через мембрану)).
Также я посмотрела лучшую по e-value находку (с лучшей по bit score она не совпадает) - это AL123456.3_49927 (DNA adenine methylase) с координатами 379226-379414. Для нее я нашла пересекающийся CDC CCP43041.1. Увы, этот CDC кодирует неизвестный белок neighborseval.tsv.
Наконец, я попробовала найти CDS моей находки с помощью поиска по аннотации кодирующих участков в геноме (по EC-кодам ферментов
- 2.1.1.37 (m5C), 2.1.1.72 (m6A), 2.1.1.113 (m4C)):
elink -target protein -db nuccore -id AL123456.3 | efilter -query '2.1.1.37' | efetch -format 'fasta'
В результате, я с помощью этой и аналогичных команд я получила 0 находок. Значит, таким методом искать не стоит.