EMBOSS, Entrez Direct, NCBI Datasets

Для выполнения практикума 9 использовали бактерию Actinomyces bowdenii

Этап 1

Идентификатор протеома в UniProt: UP000271272

Код доступа геномной сборки из GenBank: GCA_003860075.1

Код доступа геномной сборки из RefSeq: GCF_003860075.1

Ссылка на страницу сборки в Genome

Этап 2

Была скачена последовательность и feature table с помощью команды: 

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

Далее полученный файл распаковали с помощью команды: 

unzip ncbi_dataset.zip -d genome_files

Этап 3

Для корректного поиска ORF необходимо знать, какую таблицу трансляции использует исследуемый организм. Эта информация содержится в записи о таксоне в базе NCBI Taxonomy. Запись была получена с помощью EDirect следующей команды:

efetch -format xml -db taxonomy -id 131109 > taxon.xml

В файле taxon.xml найдено поле "11". Значит наш организм Actinomyces bowdenii использует 11 таблицу кода.

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

getorf -sequence ncbi_dataset/data/*/*.fna -outseq intermed.faa -table 11 -minsize 150
getorf – позволяет найти открытые рамки считывания в нуклеотидной последовательности.
intermed.faa - выходной файл, сюда getorf запишет все найденные аминокислотные последовательности ORF.
table 11 - номер таблицы генетического кода, который getorf должен использовать для перевода нуклеотидов в аминокислоты.
minsize 150 - минимальная длина ORF в нуклеотидах (50 амк), которую программа будет выводить.

Чтобы убедиться, что мы все правильно отсортировали, выполнили команду:

infoseq -only -length intermed.faa | awk '$1 < 50'
infoseq – позволяет получить базовую информацию про последовательности в виде таблицы.
-only -length — выводит только длину.
awk '$1 < 50'— первая колонка (длина последовательности) меньше 50.

Команда ничего не выдала, значит мы все сделали правильно, можем двигаться дальше.

Для создания белковой базы proteome для blastp была использована следующая команда: 

makeblastdb -in intermed.faa -dbtype prot -out proteome

Этап 4

Далее были скачены три известные ДНК-метилтрансферазы из Swiss-Prot (P0AED9, P0AEE8 и P23941), чтобы потом использовать их как запросы для поиска по сходству.

Команда: 

 echo "sw:P0AED9 sw:P0AEE8 sw:P23941" | tr " " "\n" | seqret @/dev/stdin -filter > query.fasta
seqret @/dev/stdin - программа seqret берёт свой входной список (ID белков) не из файла, а прямо из того, что ей передали через |
filter - заставляет seqret выводить результат в stdout, чтобы перенаправить его в файл через >.

Этап 5

После скачивания белков (этап 4) был запущен локальный blastp с выдачей в формате 7 (табличный формат) с помощью команды:

blastp -query query.fasta -db proteome -outfmt 7 -out blast_results.txt
  • Выдача

    • По результатам поиска во всех трёх группах были обнаружены совпадения, однако наибольший вес (bitscore) показал хит для m5C-метилтрансферазы Dcm (P0AED9). Лучшей по весу оказалась рамка NZ_RQZC01000032.1_237 (getorf указал координаты [7167–6811] по минус-цепи внутри данного контига) с bitscore = 30.4. Таким образом, данная ORF является наиболее вероятным гомологом m5C-ДНК-метилтрансферазы.

    Чтобы найти соседние гены рядом с нашей ORF, был создан файл CDS.tsv, куда были вынесены координаты всех аннотированных CDS из файла genomic.gff. Для этого выбрали только строки типа CDS и оставила нужные столбцы (4, 5, 7 и 9) с помощью команды:

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

    Для определения ближайших аннотированных CDS к найденной ORF была использована команда:

    echo -e "6811\t7167\t-\tFOUND-ORF" | cat - CDS.tsv | sort -n | grep -C 3 "FOUND-ORF" > neighbors.tsv
  • Выдача

    • В файле neighbors.tsv присутствует только сама найденная ORF, без соседних CDS. Это означает, что в диапазоне координат 6811–7167 нет аннотированных генов, и ORF не пересекается с CDS.

    Этап 6

    Для проверки, можно ли найти соответствующую МТазу по аннотациям. Был выполнен поиск белков, связанных с нуклеотидной последовательностью NZ_RQZC01000032.1, по EC-коду m5C-метилтрансферазы (2.1.1.37): 

    elink -id NZ_RQZC01000032.1 -db nuccore -target protein | efilter -query "(EC=2.1.1.37)" | efetch -format acc
    elink - получение записей по ссылкам из других записей
    db nuccore - база нуклеотидов (геномные последовательности).
    efetch -format acc - выводит только accession-номера найденных белков

    Поиск не вернул ни одной записи. Следовательно, по аннотации найти эту МТазу невозможно, так как EC-код для соответствующего гена отсутствует.