учебная страничка Маши Смирновой

практикум 9

этап 1: получение идентификаторов геномной сборки и таксона

в прошлом семестре я работала с референсным протеомом организма Halobellus clavatus, его данные:

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

загрузить нуклеотидную последовательность генома и таблицу признаков в формате gff3.

datasets download genome accession GCF_900107195.1 --include gff3,genome
unzip ncbi_dataset.zip

пояснение

  • datasets download genome accession - официальный инструмент ncbi для загрузки данных генома.
  • GCF_900107195.1 - код нашего генома.
  • --include gff3,genome - просим скачать и "карту" генома (gff3), и сам текст днк (genome).
  • unzip - распаковывает скачанный архив.

этап 3: поиск и трансляция открытых рамок считывания (orf)

сначала надо определить генетический код организма, чтобы искать открытые рамки считывания, и найти все участки днк, которые могут быть белками (не короче 50 аминокислот).

1.

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

пояснение:

  • efetch - достает (fetch) данные из базы ncbi.
  • -db 'taxonomy' - смотрим в базу систематики организмов.
  • -id '660517' - номер нашего организма в этой базе.
  • -format 'xml' - выдаем ответ в виде формата xml

результат: организм использует таблицу №11 (bacterial, archaeal and plant plastid) то есть:

<geneticcode>
    <gcid>11</gcid>
    <gcname>bacterial, archaeal and plant plastid</gcname>
</geneticcode>

как и большинство бактерий, архей и пластид

2. поиск orf программой getorf - ищем открытые рамки считывания:

getorf ncbi_dataset/data/GCF_900107195.1/*.fna -outseq orfs.orf -table 11 -minsize 150 -find 0

разбор getorf:

  • getorf - ищет открытые рамки считывания (orf).
  • -outseq - название файла с результатом.
  • -table 11 - используем найденный ранее код архей.
  • -minsize 150 - ищем рамки от 150 нуклеотидов (чтобы белок был от 150/3=50 аминокислот).
  • -find 0 - ищем открытые рамки считывания между стоп-кодонами.

проверка на то, что в каждой orf находится не менее 50 а.о.

infoseq orfs.orf -only -length | sort -n | head

пояснение

  • infoseq - извлекает информацию о длине каждого белка.
  • sort -n | head - сортирует по убыванию и показывает первые строки

проверка подтвердила, что в каждой orf не менее 50 а.о.

3. создание базы

makeblastdb -in orfs.orf -dbtype prot -out proteome

пояснение:

  • makeblastdb - превращает текстовый файл в базу данных, в которой программа blast сможет быстро искать.
  • -in orfs.orf - файл с подходящими orf
  • -dbtype prot - белковая база данных
  • -out proteome - имя бд

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

теперь нам надо скачать последовательности белков-эталонов (m5C, m6A, m4C) из swiss-prot.

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

разбор команды:

  • echo -e - печатает список идентификаторов белков.
  • sw: - показывает, что белок нужно искать в базе данных swiss-prot.
  • seqret - достает последовательности из баз.
  • @stdin - говорит программе читать список id через конвейер из команды
  • query.fasta - имя файла

этап 5: поиск по сходству последовательностей (blast)

сейчас мы сравниваем эталоны с нашими orf и найти самого похожего кандидата.

blastp -query query.fasta -db proteome -out mtf.out -outfmt 7

разбор blastp:

  • blastp - сравнивает белок с белком.
  • -query - наш эталон.
  • -db - наша база данных orf.
  • -out mtf.out - имя файла
  • -outfmt 7 - табличный формат с комментариями
mtf.out
результат:
DCM_ECOLI NZ_FNPB01000002.1_173 28.293 205 118 5 86 290 46 221 7.08e-13 68.9
  • находка: рамка NZ_FNPB01000002.1_173
  • e-value: 7.08e-13 - лучший результат
  • вес: 68.9
  • гомологом обнаружена m5C-мтаза (dcm).

координаты находки:

grep "NZ_FNPB01000002.1_173" orfs.orf
>NZ_FNPB01000002.1_173 [30779 - 31936] Halobellus clavatus strain CGMCC 1.10118, whole genome shotgun sequence

анализ окружения

с помощью команд grep, cut и sort мы выделили аннотированные гены рядом с нашей находкой.

grep "NZ_FNPB01000002.1" genomic.gff | grep -P "\tCDS\t" | cut -f4,5,7,9 > CDS.tsv

echo -e '30779\t31936\t-\tFOUND-ORF' | cat - CDS.tsv | sort -n | grep -C 3 'FOUND-ORF' > neighbors.tsv
  • -p - для понимания символа табуляции.
  • \tcds\t - мы ищем слово cds, окруженное знаками табуляции для обнаружения именно записи в третьей колонке

наша находка (30779-31936) практически совпадает по координатам с официальным геном WP_089765343.1 (30887-31939), который аннотирован как днк-цитозиновая метилтрансфераза. - dna cytosine methyltransferase (WP_089765343.1) dna (cytosine-5-)-methyltransferase activity (m5C).

neighbors.tsv

этап 6: поиск по аннотациям через e-direct

теперь надо найти метилтрансферазы в официальной аннотации ncbi, используя ферментативные коды (ec 2.1.1.37, 2.1.1.72, 2.1.1.113).

команда (для m5C и других по аналогии):

elink -db nuccore -id NZ_FNPB01000002.1 -target protein | efilter -query "2.1.1.37 [ECNO]" | efetch -format acc

пояснения:

  • elink - связывает записи в разных базах данных ncbi.
  • -db nuccore - поиск из базы nucleotide
  • -id NZ_FNPB01000002.1 — указываем конкретный контиг
  • -target protein - берет все белки закодированные на этом днк
  • efilter - фильтрует белки
  • -query "2.1.1.37 [ECNO]" оставляем только те белки, у которых в поле ecno стоит код 2.1.1.37.
  • efetch -format acc - печатает только краткий номер белка.

результат: найдено совпадение для m5C - белок WP_089765343.1. по остальным кодам (m6A и m4C) записей не найдено.