Для получения TaxID организма, которому принадлежит геном, был использован следующий конвейер:
grep 'TaxID' UP000681035.swiss | tr ' ' '\n' | grep 'TaxID' | tr '=' '\n' | grep -v 'TaxID' | sort -u
Результат: 2714355
Команда для получения списка сборок, соответствующих таксону:
datasets summary genome taxon 2714355 --as-json-lines | dataformat tsv genome
Я получила несколько сборок и, чтобы найти среди них ту, которая соотвествует протеому моей бактерии, воспользовалась конвейером:
zgrep 'DR EMBL' UP000681035.swiss | cut -f 5 -d ' ' | tr -d ';' | elink -db nuccore -target assembly | esummary | xtract -pattern DocumentSummary -element AssemblyAccession
Нужная сборка: GCF_018408575.1
Я загрузила последовательность и feature table с помощью команды:
datasets download genome accession GCF_018408575.1 --include gff3,genome
Чтобы определить, какой вариант генетического кода использует исследуемый мной организм, я ввела конвейер:
efetch -db 'taxonomy' -id '2714355' -format xml
В полученном файле с выдачей написано, что бактерия использует 11-ую таблицу генетического кода. Далее я нашла открытые рамки считывания:
getorf GCF_018408575.1_ASM1840857v1_genomic.fna ORFs.fasta -minsize 150 -table 11 -filter
Чтобы проверить длину наименьшей из полученных последовательностей, я использовала конвейер:
infoseq -only -length ORFs.fasta | sort -n | head -n 2
Далее я создала белковую базу данных на основе полученных открытых рамок считывания:
makeblastdb -in ORFs.fasta -out ORFs -dbtype prot
Я скачала последовательности белков: P0AED9 (Dcm, m5C-МТаза, E.coli), P0AEE8 (Dam, m6A-МТаза, E.coli), и P23941 (m4C-МТаза, Bacillus amyloliquefaciens):
echo -e 'sw:P0AED9\nsw:P0AEE8\nsw:P23941' > s.txt | seqret @s.txt query_MTases.fasta
С помощью blastp на основе созданной базы данных ORFs и скачанных последовательностей ДНК-метилтрансфераз я осуществила поиск схожих фрагментов белковых последовательностей:
blastp -query query_MTases.fasta -db ORFs -out blastp.out -outfmt 7
Лучшая по весу находка - NZ_AP023418.1_4753, ее координаты - 950508 - 951416, вес - 123, она является гомологом DNA adenine methylase (m6A).
Команда для создания файла CDS.tsv, содержащего информацию о CDS, находящихся в той же геномной последовательности, что и интересующая находка:
grep 'NZ_AP023418.1' genomic.gff | grep CDS | cut -f4,5,7,9 > CDS.tsv
Чтобы получить находки с близкими к лучшей находке координатами, я ввела конвейер:
echo -e '950508\t951416\t-\tMY-ORF' | cat - CDS.tsv | sort -n | grep -C 3 'MY-ORF' > neighbors.tsv
У меня получилось, что один CDS перекрывается с лучшей находкой. Его координаты - 950520 - 951419, они не совпадают с координатами лучшей находки.
Я попробовала найти CDS, соответствующий моей находке, с помощью поиска по аннотации кодирующих участков в геноме. Для этого я написала команду:
elink -target protein -db nuccore -id NZ_AP023418.1 | efilter -query '2.1.1.72' | efetch -format 'fasta'
Был найден 1 белок, тот же, который я нашла в предыдущем пункте - ДНК-метилтрансфераза m6A.