Задание 0
В прошлом семестре мной был выбран протеом UP000239526 бактерии Providencia rettgeri. Файл с белками протеома в формате swiss у меня сохранился, поэтому скачивать заново его не пришлось.
Задание 1. Получение AC геномной сборки
Найдём TaxID организма:
grep '^OX' UP000239526.swiss| cut -c -19 | sort | uniq -c
Выдача:
4034 OX NCBI_TaxID=587
Посмотрим на возможные геномные сборки для данного организма:
datasets summary genome taxon 587 --as-json-lines | dataformat tsv genome | less
AC: GCF_013702265.1
Задание 2. Cкачивание последовательности генома и таблицы локальных особенностей
При помощи следующего алгоритма скачаем геном и таблицу локальных особенностей (опция include gff3,genome добавляет таблицу локальных особенностей и последовательность генома в архив):
datasets download genome accession GCF_013702265.1 --include gff3,genome
Команда для распаковки архива:
unzip ncbi_dataset.zip
Задание 3. Поиск и трансляция открытых рамок считывания
Для начала определим какой вариант генетического кода использует Providencia rettgeri:
efetch -db 'taxonomy' -id '587' -format 'xml'|less
Находим строчку, начинающуюся с GCId: < GCId>11< /GCId>, следовательно, в моём случае используется таблица генетического кода №11.
Найдём открытые рамки считывания (-table 11 - задает таблицу №11 для трансляции; -minsize 150 - минимальный размер открытых рамок):
getorf -sequence 'GCF_013702265.1_ASM1370226v1_genomic.fna' -minsize 150 -table 11 out.fasta
Где out.fasta - файл с искомыми открытыми рамками считывания, а GCF...fna - нуклеотидная последовательность генома. Затем проверим, что среди транслированных белков нет элементов короче 50 а.о., с помощью алгоритма:
infoseq out.fasta -only -length | sort -n | head -2 | less
Элементов длиной меньше 50 не найдено, поэтому дальше создадим белковую базу данных на основе полученных открытых рамок считывания:
makeblastdb -dbtype prot -in out.fasta -out ORFs
Задание 4. Получение последовательностей гомологичных метилтрансфераз
Скачаем последовательности белков разных ДНК-метилтрансфераз: P0AED9 (Dcm, m5C-МТаза, E.coli), P0AEE8 (Dam, m6A-МТаза, E.coli) и P23941 (m4C-МТаза, Bacillus amyloliquefaciens).
echo -e 'sw:P0AED9\nsw:P0AEE8\nsw:P23941'|seqret -filter 'list::stdin' -outseq 'query_MTases.fasta'
Задание 5. Поиск по сходству последовательностей
Полученные последовательности метилтрасфераз используем в роли запроса для blatp по созданной ранее базе данных ORFs:
blastp -query query_MTases.fasta -db ORFs -outfmt 7 > out.txt
Выдача программы blastp находится в файле. Самая лучшая находка получилась с с m6A-МТазой из E.coli, её вес - 389; название рамки - NZ_CP059348.1_22687, её координаты в геноме: 129349 - 130188, E-value = 2.65e-138.
У этой находки действительно хорошие показатели, поэтому будем считать, что этот участок действительно может быть гомологичен.
Команда для нахождения координат:
grep 'NZ_CP059348.1_22687' out.fasta
Теперь при помощи координат мы можем определить с какими кодирующими белки участками (CDS) пересекается моя лучшая находка, для этого создадим вспомогательный файл CDS.tsv, содержащий только CDS из той же нуклеотидной последовательности, что и моя находка, а также их координаты и аннотацию, все данные о CDS возьмём из ранее скачанной таблицы локальных особенностей (genomic.gff):
grep 'NZ_CP059348.1' genomic.gff | grep CDS | cut -f4,5,7,9 > CDS.tsv
Найдём участки близкие/перекрывающиеся с моей находкой при помощи конвейера:
echo -e '129349\t130188\t-\tMY-ORF' | cat - CDS.tsv | sort -n | grep -C 3 'MY-ORF' > neighbors.tsv
Нашлась CDS - аденинспецифичная ДНК-метилтрансфераза (ID=cds-WP_154602134.1). Координаты: 129346 - 130182, CDS полностью лежит внутри моей находки.
Задание 6. Поиск по аннотациям кодирующих участков
Проверим, возможно ли найти нужную нам CDS по аннотациям в геноме.
Запустим конвейер для поиска с использованием кода класса найденного фермента:
elink -target protein -db nuccore -id NZ_CP059348.1 | efilter -query '2.1.1.72' | efetch -format 'fasta' | less
Получилось 5 находок, в числе которых оказалась найденная мной в предыдущем пункте аденинспецифичная ДНК-метилтрансфераза (WP_154602134.1).
Идентификаторы остальных найденных белков: WP_025611180.1, WP_166266817.1, WP_181478041.1, WP_181478211.1 (тоже метилтрансфераза).