Практикум 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) записей не найдено.