На первом курсе моим объектом была Pussilibacter faecalis, в том числе для анализа протеома. В данной работе нам понадобятся следующие данные
Для того, чтобы скачать последовательность и таблицу нужен код доступа RefSeq и NCBI Datasets command-line tool. Вводим команду:
datasets download genome accession GCF_018408705.1 --include genome,gff3 --filename genome_data_pus.zip
Затем с помощью
unzip genome_data_pus.zip
Распаковываем zip файл, полученный выше.
С помощью efetch можно скачать информацию про наш таксон в формате xml, чтобы посмотреть, какой вариант генетического кода использует организм.
efetch -db taxonomy -id 2714358 -mode xml > pusib.xml
У нас используется таблица 11 (бактериальная). По файлу генома ищем рамки между двумя стоп-кодонами, между которыми не менее 50 аминокислот. С помощью этой команды находим такие рамки и транслируем их:
getorf -table 11 -minsize 150 -find 0 ncbi_dataset/data/GCF_018408705.1/GCF_018408705.1_ASM1840870v1_genomic.fna puspro.fa
Проверим нет ли среди последовательностей тех, которые короче 50 аминокислот:
infoseq -only -length puspro.fa | sort | uniq
Среди чисел нет тех, которые меньше 50. Также подготовим на основе этого файла протеом для blastp:
makeblastdb -in puspro.fa -dbtype prot -out proteome
Нужно получить последовательности P0AED9 (Dcm, m5C-МТаза, E.coli), P0AEE8 (Dam, m6A-МТаза, E.coli), и P23941 (m4C-МТаза, Bacillus amyloliquefaciens). Это можно сделать с помощью seqret:
echo "sw:P0AED9 sw:P0AEE8 sw:P23941" | tr ' ' '\n' | seqret @stdin -outseq query.fasta
Дальше ищем с помощью blastp гомологичные последовательности среди ранее полученного протеома для трех метилтрансфураз, которые мы тоже нашли ранее:
Получаем файл в виде таблицы. Самая лучшая по весу (67.0) находка это NZ_AP023421.1_47 для Dam, m6A-МТаза, E.coli. Эта последовательность кодируется в pMM59_01 плазмиде с координатами [9185 - 11050] на + цепи.
Затем из таблицы локальных особенностей находим, какие CDS располагаются рядом. Для этого отберем сначала все кодирующие последовательности из плазмиды, их координаты, цепью и доп. информацией:
Затем как раз определим те, которые находятся рядом:
echo -e '9200\t11053\t-\tFOUND-ORF' | cat - CDS.tsv | sort -n | grep -C 3 'FOUND-ORF' > neighbors.tsv
Для поиску по ЕС ферментов, я ввел следующие команды:
Не было найдено ни одной находки.