Для выполнения этого задания был выбран протеом Brachybacterium nesterenkovii (UP000195981). Сначала я получила TaxID организма, которому принадлежит геном, использовав конвейер:
grep "^OX" UP 000195981.swiss | tr ' ' '\n'|grep "^NCBI_TaxID=" | sort -u
Я получила TaxID организма: NCBI_TaxID=47847.
Далее я нашли АС сборки, получив соответствующие записи в Assembly, для этого я использовала конвейер:
grep 'DR' UP000195981.swiss | grep 'EMBL;' | cut -d ' ' -f 5 | tr -d ';'| elink -db 'nuccore' -target 'assembly ' | esummary | xtract -pattern DocumentSummary -element AssemblyAccession
Полученный результат содержал АС: GCF_900163655.1
Имея AC геномной сборки, я загрузила геном и таблицу локальных особенностей с помощью:
datasets download genome accession GCF_900163655.1 --include gff3 --include genome
Параметр --include gff3 указывает на добавление к архиву таблицы локальных особенностей в формате GFF3
После загрузки я распаковала архив с файлами: unzip ncbi_dataset.zip
Для поиска открытых рамок считывания я определила вариант генетического кода, который использует исследуемый организм:
efetch -db 'taxonomy' -id '47847' -format 'xml' | less
Элемент
Далее я нашла открытые рамки считывания, использовав конвейер:
getorf -sequence 'GCF_900163655.1_ASM90016365v1_genomic.fna' -outseq 'out' -table 11 -minsize 150
Параметр -table указывает на необходимую таблицу генетического кода, а параметр -minsize указывает на минимальный размер открытых рамок.
По полученным открытым рамкам я создала по ним белковую базу с названием 'ORFs': makeblastdb -in 'out' -dbtype prot -out 'ORFs'
Я проверила отсутствие элементов в трансляциях короче 50 аминокислотных остатков с помощью:
infoseq 'out' -only -length|sort -n|head -2|less
Из базы данных Swiss-Prot я скачала три последовательности различных ДНК-метилтрансфераз (P0AED9 (Dcm, m5C-метилтрансфераза, Escherichia coli), P0AEE8 (Dam, m6A-метилтрансфераза, Escherichia coli) и P23941 (m4C-метилтрансфераза, Bacillus amyloliquefaciens)), сгенерировав один файл (query_MTases.fasta) и использовав конвейер:
echo -e 'sw:P0AED9\nsw:P0AEE8\nsw:P23941'|seqret -filter 'list::stdin' -outseq 'query_MTases.fasta'
Я провела поиск по сходству последовательностей:
blastp -query 'query_MTases.fasta' -db 'ORFs' -outfmt 7 -out 'output_table'
По результатам поиска я получила файл 'output_table' со следующей информацией:
Находок с DCM_ECOLI: 1
Нахдок с DMA_ECOLI: 3
Находок с MTB1_BACAM: 0
Проанализировав находки, лучшей была выбрана находка с m5C-МТаза из E.coli (NZ_FWFG01000060.1_90 – название рамки), так как ее вес составил 27.7, наибольший из всех находок, и E-value равно 6.8. У этой находки высокое значение E-value и малый вес, что позволяет предположить, что этот участок может быть не гомологичным. Его координаты в геноме: 17762 - 17911 (данные о координатах рамки содержатся в файле 'out')
Я определила, какие CDS из таблицы локальных особенностей генома располагаются рядом, выбрав необходимые столбцы и используя конвейер:
grep 'CDS' genomic.gff | grep '^NZ_FWFG01000060.1' | cut -f 4,5,7,9>CDS.tsv
Далее я выбрала соседние к выбранной находке кодирующие последовательности:
echo -e '17762\t17911\t-\tMY-ORF' | cat - CDS.tsv | sort -n | grep -C 3 'MY-ORF' > neighbors.tsv
В полученном файле была найдена CDS (ID=cds-WP_087103857.1), которая является carbohydrate ABC transporter permease. Координаты находки (17762 - 17911) и CDS (16878 - 17858) пересекаются.
Я проверила возможность нахождения CDS соответствующего находке с помощью поиска по аннотации кодирующих участков в геноме, используя конвейер:
elink -target 'protein' -db 'nuccore' -id ''NZ_FWFG01000060.1' | efilter -query '2.1.1.37' | efetch -format 'fasta' | less
В результате ни одного белка по данному запросу не было найдено.