Работа с NCBI при помощи Entrez Direct.
И снова EMBOSS.
Получение AC геномной сборки и TaxID организма.
Для штамма ATCC 19860 бактерии Paracidovorax avenae (ASM17685v2) были получены идентификаторы в базах данных:
Uniprot proteomes: UP000002482
Геномная сборка GenBank: GCA_000176855
Геномная сборка RefSeq: GCF_000176855.2
Скачивание данных с NCBI
Для последующих заданий была скачена последовательность единственной хромосомы бактерии и таблицу локальных особенностей в формате GFF3.
datasets download genome accession GCF_000176855.2 --include gff3,genome
Поиск и трансляция открытых рамок считывания
При помощи средств EMBOSS нам нужно получить все открытые рамки считывания по скаченной в предыдущем пункте последовательности с условием, что мы анализируем рамки считывания длинной от 150 нуклеотидов. Но для того, чтобы это сделать правильно, нужно знать, какая таблица генетического кода используется этой бактерией. Это информация сожержится в записи генома в базе данных Nucleotide на NCBI:
elink -db assembly -id GCF_000176855.2 -target nuccore | esummary | xtract -pattern 'DocumentSummary' -element 'GeneticCode' | uniq > GeneticCodeTable.txt
Файл на выходе содержал число 11 - такой номер имеет таблица генетического кода имеет данный организм. Теперь при помощи программы getorf из пакета EMBOSS мы можем получить последовательности всех белков, с распознанных рамок считывания:
getorf -sequence GCF_000176855.2_ASM17685v2_genomic.fna -table 11 -minsize 150 -outseq proteins
Стоит проверить, нигде ли мы не ошиблись - правильную ли минимальную длину мы указали. Так как минимальная длина рамки считывания была равна 150, следовательно минимальная длина белка должна быть больше либо равна 50 ак. Проверим это при помощи программы infoseq из пакета EMBOSS:
infoseq -sequence proteins -only -length -noheading | sort -n | head -n 1 | less
Программа выводит 50, что означает, что все выполнено верно.
Теперь на основе полученного списка белков нужно сделать базу данных для дальнейшей работы локальной программы Blast, а именно blastp:
makeblastdb -in proteins -dbtype prot -out proteome
Получение последовательностей гомологичных метилтрансфераз
В качестве организма из которого мы берем последовательности метилтрансфераз, по которым мы будем находить гомологичные в Paracidovorax avenae была выбрана E.coli. При помощи программы seqret из пакте EMBOSS последовательности были скачаны в файл query.fasta:
echo -e "sw:P0AED9\nsw:P0AEE8\nsw:P23941" | seqret -sequence @stdin -outseq query.fasta
Поиск по сходству последовательностей
При помощи локального blastp мы искали схожие с метилтрансферазой последовательности среди белков Paracidovorax avenae в ранее полученной базе данных proteome. Файл был получен при помощи следующей команды:
blastp -query query.fasta -db proteome -out results.txt -outfmt 7 -max_target_seqs 1
Выходной файл: results.txt. Как можно заметить, в команде использовалась опция -max_target_seqs, что позволяет показывать лишь первую находку среди всех в списке. В выходном файле содержатся выводы для всех трех последовательностей метилтрансфераз E.Coli.
Лучшей находкой является находка при поиске по последовательности ДНК-цитозин метилтрансферазы(2.1.1.37) (e-value = 1.44e-11, e-value у остальных находок >0.2).
Дальше нам нужно было понять, какие аннотированные гены находятся рядом с нашей находкой.
При помощи следующей команды мы получили координаты лучшей находки:
<proteins grep 'NC_015138.1_27371' | less
Далее нам нужно было получить таблицу с CDS анализируемой бактерии:
<genomic.gff grep 'CDS' | cut -f 4,5,7,9 >CDS.tsv
Имея координаты (5065581 - 5063977) и таблицу CDS, мы осуществили поиск CDS, соседствующих с координатами нашей находки:
echo -e '5063977\t5065581\t-\tFOUND-ORF' | cat - CDS.tsv | sort -n | grep -C 3 'FOUND-ORF' > neighbors.tsv
Выходной файл: neighbors.tsv.
При анализе таблицы было замечено, что сильнее остальных пересекается ген с координатами 5063974-5065458, почти полностью попадающий в последовательность находки. Это ген белка ДНК-цитозин метилтрансферазы. Ее мы, собственно, и искали.
Поиск по аннотациям кодирующих участков
В заключение, интересно узнать, смогли ли бы мы найти этот белок по аннотированным генам бактерии Paracidovorax avenae. Для этого был создан запрос-конвейер entrez:
elink -db nuccore -id 'NC_015138.1' -target 'protein' | efilter -query '2.1.1.37[EC/RN Number]' | efetch -format 'ipg' | less
На выходе мы получаем запись в базе данных NCBI Protein. ID белковой записи - WP_244875500.1. А координаты, указанные в этом формате, совпадают с координатами из предыдущего пункта (5063974-5065458).
То есть в данном случае, мы могли бы найти запись метилтрансферазы просто при помощи поиска по аннотацияем. Однако, можно сказать, что нам повезло, так как записи могло бы и не быть, а первый способ поиска генов белков со схожими доменами никак не зависит от ее наличия.