Практикум 9. Поиск ДНК-метилтрансфераз в геноме археи по последовательности и по аннотации

Этап 0: протеом археи

Белки протеома археи Natronorubrum halalkaliphilum были получены в практикуме 8 второго семестра следующей командой:

wget 'https://rest.uniprot.org/uniprotkb/stream?compressed=true&format=txt&query=proteome:UP000434101' -O UP000434101.swiss.gz

Этап 1: получение AC геномной сборки

Через поиск в UniProt Proteoms и NCBI Assembly (Datasets Genome) была получена следующая информация.

Идентификатор протеома Proteome ID - UP000434101
Ссылка на запись в Proteoms: https://www.uniprot.org/proteomes/UP000434101

RefSeq идентификатор Natronorubrum halalkaliphilum - GCF_009834785.1

Соответствующий код доступа геномной сборки из GenBank - GCA_009834785.1
Ссылка на страницу сборки в Genome: https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_009834785.1/

Этап 2: скачивание последовательности генома и таблицы локальных особенностей

Последовательность и feature table были загружены с помощью команды:

datasets download genome accession GCF_009834785.1 --include genome,gff3

И далее распакованы.

Этап 3: поиск и трансляция открытых рамок считывания

Прежде чем искать открытые рамки считывания, нужно определить, какой вариант генетического кода использует исследуемый организм. Эту информацию получили из записи про таксон в базе NCBI Taxonomy, скачав её с помощью efetch в формате xml:

efetch -db 'taxonomy' -id '2691917' -format 'xml'

В поле GeneticCode указано, что наша архея использует таблицу 11.

Найти открытые рамки считывания и сразу получить их трансляции можно с помощью программы getorf:

getorf -sequence 'GCF_009834785.1_ASM983478v1_genomic.fna' -table 11 -minsize 150 -find 0 'database.fasta'

Получаем рамки между стоп-кодонами (чтобы не потерять фрагмент последовательности из-за неверного старт-кодона), при этом их трансляции не короче 50 аминокислот. Последовательности трансляций сохранили в промежуточный файл. Длину проверили следующей командой:

infoseq 'database.fasta' -only -length | tail -n +2 | sort -nu | less

Создали по ним белковую базу для blastp, назвали её proteome:

makeblastdb -in 'database.fasta' -dbtype 'prot' -out 'proteome'

Этап 4: получение последовательностей гомологичных метилтрансфераз

Получаем последовательности следующих белков (указаны коды доступа в базе Swiss-Prot): P0AED9 (Dcm, m5C-МТаза, E.coli), P0AEE8 (Dam, m6A-МТаза, E.coli), и P23941 (m4C-МТаза, Bacillus amyloliquefaciens).

echo 'sw:P0AED9 sw:P0AEE8 sw:P23941' | tr ' ' '\n' | seqret @stdin -filter 'query.fasta'

Этап 5: поиск по сходству последовательности

С помощью локального blastp в нашей локальной базе proteome и query.fasta в качестве запроса были найдены гомологи метилтрансфераз из этапа 4.

blastp -db 'proteome' -query 'query.fasta' -outfmt 7 -out 'MTas.tab'

Полученный файл: MTas.tab

Лучшей по весу находкой (score = 54.7) оказалась гомологичная m6A-МТазе E.coli Deoxyadenosyl-methyltransferase с названием рамки NZ_WUYX01000053.1_2398. Координаты рамки в геноме: [43634- 44569], - были определены следующей командой:

grep -e 'NZ_WUYX01000053.1_2398' database.fasta

Маленький вес находки обусловлен тем, что мы сравниваем архею с бактерией.

Ищем, какие CDS из таблицы локальных особенностей генома располагаются рядом.

grep 'NZ_WUYX01000053.1' genomic.gff | cut -f 3,4,5,7,9 | grep '^CDS' | cut -f 2-5 > CDS.tsv

Далее используем следующий конвейер для поиска CDS из таблицы локальных особенностей генома, располагающихся рядом с нашей находкой:

echo -e '44569\t43634\t-\tFOUND-ORF' | cat - CDS.tsv | sort -n | grep -C 3 'FOUND-ORF' > neighbors.tsv

Одна из соседних CDS пересекается с нашей рамкой, она имеет координаты [43631 - 44557]. В описании данной CDS указано site-specific DNA-methyltransferase (adenine-specific) activity, - что подтверждает гомологичность нашей и найденной последовательности(идентификатор - WP_160066311.1).

Ссылка на файл: neighbors.tsv

Этап 6: поиск по аннотациям кодирующих участков

В завершении, определяем CDS, соответствующий нашей находке с помощью поиска по аннотации кодирующих участков в геноме. Поиск проводим по аннотации и по EC-коду m6A ДНК-метилтрансферазы в базе nuccore.

elink -db nuccore -id 'NZ_WUYX01000053.1' -target protein | efilter -query '2.1.1.72' | efetch -format 'acc'

Elink ищет связанные записи между базами данных. Он отбирает нашу нуклеотидную последовательность в базе Nucleotides по ее id и находит по ссылкам все связанные с ней белковые записи из базы Protein. Далее efilter фильтрует по нужному нам ЕС-коду m6A-МТазы(список полей можно получить с помощью einfo -db 'protein' -fields). Следом выводятся только Accession коды белковых последовательностей.

В выдаче получаем одну запись WP_160066311.1, чей идентификатор совпадает с нашей находкой из предыдущего пункта.