Практикум 9. EMBOSS, Entrez Direct, NCBI Datasets

Протеом бактерии Thiomicrorhabdus sediminis (G1) (ID штамма: 2580412) сохранился с прошлого семестра (Практикум 8). В дальнейшем я буду работать с ним (и списывать у самого себя из прошлого).

Этап 1. Получение AC геномной сборки и TaxID организма

Таблица 1. Информация о геномных сборках и TaxID организма

ID сборки на RefSeq	GCF_005885815.1
ID сборки на INSDC	GCA_005885815.1
ID протеома на UniProt Proteomes	UP000304864
TaxID	2580412

Были представлены геномные сборки как в GenBank, так и в RefSeq. В дальнейшем будет использоваться сборка из RefSeq (GCF_005885815.1).

Этап 2. Скачивание последовательности генома и таблицы локальных особенностей

Использованные команды:
datasets download genome accession GCF_005885815.1 --include gff3,genome
unzip ncbi_dataset.zip
Для удобства работа велась в папке public_html/term3/pr9/dataset/.

Этап 3. Поиск и трансляция открытых рамок считывания

С помощью команд
efetch -db 'taxonomy' -id '2580412' -format 'xml' > entry.xml
cat entry.xml | grep 'GCId'
было получено, что Thiomicrorhabdus sediminis использует 11 таблицу генетического кода.

Поиск и трансляция открытых рамок считывания (Open Reading Frame, ORF) были выполнены с помощью команды
getorf GCF_005885815.1_ASM588581v1_genomic.fna -outseq translated_orfs.fasta -table 11 -minsize 150 -find 0
Объяснение команды:

GCF_005885815.1_ASM588581v1_genomic.fnaa – fasta-файл с нуклеотидной последовательностью генома;
translated_orfs.fasta – файл с аминокислотной последовательностью транслированных ORF;
-table 11 – указание таблицы генетического кода;
-minsize 150 – указание минимальной длины ORF в нуклеотидах (для получения последовательностей длиной не менее 50 аминокислот);
-find 0 – Значение 0 означает транслирование последовательности между стоп-кодонами (есть 2 определения ORF: область без стоп-кодонов или область между старт- и стоп-кодонами)

С помощью команды
infoseq translated_orfs.fasta -only -length | sort -n | less
было подтверждено, что транслированные ORF имеют длину не менее 50 аминокислот.

Белковая база для будущего blastp была создана с помощью команды
makeblastdb -in translated_orfs.fasta -out proteome -dbtype 'prot'

Этап 4. Получение последовательностей гомологичных метилтрансфераз

Я создал файл query.fasta, содержащий последовательности трех ДНК-метилтрансфераз: P0AED9 (Dcm, m5C-МТаза, E.coli), P0AEE8 (Dam, m6A-МТаза, E.coli), и P23941 (m4C-МТаза, Bacillus amyloliquefaciens). Использованная команда:
echo "sw:P0AED9 sw:P0AEE8 sw:P23941" | tr ' ' '\n' | seqret -filter '@stdin' -outseq query.fasta

Этап 5. Поиск по сходству последовательностей

С помощью метода blastp по локальной базе данных proteome был проведен поиск последовательностей, обладающих наибольшим сходством.
blastp -query query.fasta -outfmt 7 -db proteome -out blastp_out.txt
Файл с выдачей blastp: blastp_out.txt

Рамка NZ_CP040602.1_1456 имела наибольший вес – 75.9 (evalue: 9.43e-16). Обнаруженный гомолог: Dcm, m5C-МТаза, E.coli. Координаты: 460215 - 461321. Команда, использованная для поиска координат:
grep 'NZ_CP040602.1_1456 ' translated_orfs.fasta
Пробел после запроа в grep использовался для отсеивания результатов вида *_14560, *_14561 и т.п.

Из таблицы локальных особенностей были отобраны строки, соотвествующие CDS, и столбцы с координатами (4,5), цепью (7) и дополнительной информацией (9)
grep 'CDS' genomic.gff | cut -f1,3,4,5,7,9 > CDS.tsv
С помощью команды
cut -f2 CDS.tsv| sort -u | less
я проверил, что были отобраны только CDS длиной не менее 50 нуклеотидов. Затем я попробовал отобрать соседние с искомым CDS с помощью следующей команды:

echo -e 'NZ_CP040602.1_1456\tCDS\t460215\t461321\t-\tFOUND-ORF' | cat - CDS.tsv | sort -n | grep -C 3 'FOUND-ORF' > neighbors.tsv

Однако результат конвейера оказался неудовлетворительным (файл neighbors_old.tsv). Поэтому для поиска последовательностей, соседних с искомым CDS, я написал python-скрипт neighbor_finder.py, который сортирет исходную CSV-таблицу в порядке возрастания расстояния до заданной последовтельости (включая саму последовательность).

Мануал доступен при вызове опции -h. Команда для запуска: python .\neighbor_finder.py -start 460215 -end 461321 CDS.tsv. Файл выдачи: out.tsv. Первые пять строк выдачи:

NZ_CP040602.1	CDS	460221	461324	+	ID=cds-WP_138563936.1;Parent=gene-FE785_RS02085;Dbxref=GenBank:WP_138563936.1;Name=WP_138563936.1;Ontology_term=GO:0006304,GO:0003677,GO:0003886,GO:0008168,GO:1904047;gbkey=CDS;go_function=DNA binding|0003677||IEA,DNA (cytosine-5-)-methyltransferase activity|0003886||IEA,methyltransferase activity|0008168||IEA,S-adenosyl-L-methionine binding|1904047||IEA;go_process=DNA modification|0006304||IEA;inference=COORDINATES: protein motif:HMM:NF012372.6;locus_tag=FE785_RS02085;product=DNA cytosine methyltransferase;protein_id=WP_138563936.1;transl_table=11
NZ_CP040602.1	CDS	461345	462298	+	ID=cds-WP_138563938.1;Parent=gene-FE785_RS02090;Dbxref=GenBank:WP_138563938.1;Name=WP_138563938.1;Ontology_term=GO:0003676,GO:0004519,GO:0008270;gbkey=CDS;go_function=nucleic acid binding|0003676||IEA,endonuclease activity|0004519||IEA,zinc ion binding|0008270||IEA;inference=COORDINATES: protein motif:HMM:NF013964.6;locus_tag=FE785_RS02090;product=HNH endonuclease;protein_id=WP_138563938.1;transl_table=11
NZ_CP040602.1	CDS	458889	460121	+	ID=cds-WP_138563934.1;Parent=gene-FE785_RS02080;Dbxref=GenBank:WP_138563934.1;Name=WP_138563934.1;Ontology_term=GO:0006310,GO:0015074,GO:0003677;gbkey=CDS;go_function=DNA binding|0003677||IEA;go_process=DNA recombination|0006310||IEA,DNA integration|0015074||IEA;inference=COORDINATES: protein motif:HMM:NF012798.6;locus_tag=FE785_RS02080;product=tyrosine-type recombinase/integrase;protein_id=WP_138563934.1;transl_table=11
NZ_CP040602.1	CDS	462327	462461	+	ID=cds-WP_202978313.1;Parent=gene-FE785_RS11005;Dbxref=GenBank:WP_202978313.1;Name=WP_202978313.1;gbkey=CDS;inference=COORDINATES: protein motif:HMM:NF015790.6;locus_tag=FE785_RS11005;product=hypothetical protein;protein_id=WP_202978313.1;transl_table=11
NZ_CP040602.1	CDS	462611	462709	+	ID=cds-WP_420856750.1;Parent=gene-FE785_RS11010;Dbxref=GenBank:WP_420856750.1;Name=WP_420856750.1;gbkey=CDS;inference=COORDINATES: protein motif:HMM:NF016370.6;locus_tag=FE785_RS11010;product=hypothetical protein;protein_id=WP_420856750.1;transl_table=11

Видно, что WP_138563936.1 (460221 461324) пересекается с исследуемой CDS, однако координаты не совпадают полностью. Стоит отметить, что ее функция сходится с ожидаемой (product=DNA cytosine methyltransferase; IEA,DNA (cytosine-5-)-methyltransferase activity), что подтверждает находку

Примечания:

Программа также выводит в консоль топ 5 самых близких по расстоянию последовательностей с указанием их коордиат и расстояних до них. Параметр num – номер строки во входном файле
Написание скрипта для такой подзадачи может показаться нерациональным, однако в моем контексте это было уместно, так как мы как раз проходили модуль argparse на курсе Python, и я увидел возможность опробовать argparse в действии.
На этот раз программа была написана без использования языковых моделей, и я могу объяснить каждый элемент в скрипте

Этап 6. Поиск по аннотациям кодирующих участков

Я выполнил поиск по базе данных nuccore с помощью слеующего конвейера:

elink -db nuccore -id 'NZ_CP040602.1' -target protein| efilter -query '2.1.1.37' | efetch -format acc

Пояснение:

elink получает идентификаторы белков (-target protein) из выбранной базы данных (nuccore), которая коидируется нуклеотидной последовательностью (исходная база данных, для которой указан id)
efilter фильтруетт выдачу, используя ЕС-код фермента
efetch скачивает результат в заданном формате

Был найден белок WP_138563936.1, что означает совпадение результата с предыдущим.