Учебный сайт Левина Ильи, 3-й семестр

Команды EMBOSS, BLAST+

Задание 1

Задача 1

Проверить, сколько находок с E-value < 0.1 в среднем находит blastn для случайной последовательности данной длины в данном геноме бактерии.

Скрипт был написан на языке Python с использованием модуля subprocess, с помощью которого я мог обратиться к программам blastn, makeblastdb и makenucseq.

Здесь лежат: скрипт, fasta-файл с последовательностью генома E. coli.

Скрипт работает таким образом:

На вход подаётся 2 аргумента командной строки: первым - имя файла с геномом бактерии, вторым же нужно указать длину случайных последовательностей, которые мы будем по этому геному "бластовать";
Скрипт индексирует файл с геномом бактерии, делает из него базу для blastn;
Далее скрипт создаёт fasta-файл со 100 случайными последовательностями введённой длины;
Скрипт запускает blastn, подавая как query получившийся в пункте 3 fasta-файл, а как базу данных - проиндексированный файл с геномом бактерии;
С помощью полученного выхода blastn производится подсчёт среднего количества находок с E-value < 0.1;
Пункты 3, 4 и 5 повторяются 100 раз, и от всех полученных средних значений берётся среднее значение. Зачем? Я посчитал нужным сделать это для того, чтобы учесть ту самую "случайность" создаваемых случайных нуклеотидных последовательностей, после того, как заметил, что последовательные запуски скрипта с одними и теми же параметрами могут выдавать результаты, отличающиеся на порядок.

Пример запуска скрипта:

$ python3 script1.py Ecoli_compl_genome.fasta 25
В среднем blastn находит 0.14 находок c E-value < 0.1 для случайной последовательности длины 25 в Вашем геноме бактерии.

Таблица 1. Среднее количество находок blastn по геному E. coli при подаче на вход случайных последовательностей разной длины
Длина случ. последовательности	10	25	50	100	500	1000	5000	10000
Среднее кол-во находок	0	0.14	0.04	0.11	0.06	0.12	0.05	0.11

Как можно судить по полученным результатам, случайной последовательности может быть, конечно, найден гомолог в геноме E. coli, но это довольно редкое событие.

Задание 2

C помощью BLAST+ сделать вывод о наличии гомолога белка в неаннотированной сборке генома Amoeboaphelidium protococcarum (примитивного родственника грибов).

Для поиска гомологов я решил выбрать несколько белков, которые, я считаю, должны быть у всех эукариот:

P13382 - каталитическая альфа-субъединица ДНК-полимеразы из Saccharomyces cerevisiae. Предварительный запрос в UniProt выглядит вот так:
```
dna polymerase organism:saccharomyces cerevisiae AND reviewed:yes
```
P53942 - А-субъединица РНКазы H2 из Saccharomyces cerevisiae. Предварительный запрос в UniProt выглядит вот так:
```
ribonuclease h organism:saccharomyces cerevisiae AND reviewed:yes
```
Q04437 - 2-я субъединица АТФ-зависимой ДНК-хеликазы II из Saccharomyces cerevisiae. Предварительный запрос в UniProt выглядит вот так:
```
dna helicase organism:saccharomyces cerevisiae AND reviewed:yes
```

Далее мною был написан скрипт, который, собственно, выполняет всю техническую часть задания. Работает он так:

На вход принимается 4 аргумента: первый - имя файла с геномом организма, в котором [в геноме] мы будем искать гомологов нашим белкам; второй, третий и четвёртый - АС белков, гомологи которым мы будем искать в данном геноме;
Скрипт индексирует файл с геномом как базу данных для того, чтобы подать её BLAST'у;
Далее скрипт скачивает последовательность белка по введённому АС с помощью seqret. Шаблон команды такой:
```
$ seqret 'sw:<AC>' '<AC>.fasta' -filter
```
После этого скрипт запускает tblastn с такими параметрами:
```
$ tblastn -word_size 3 -query '<AC>.fasta' -db genome_db -max_target_seqs 10 -outfmt 7 -out 'output?.tsv'
```
где вместо "?" будут порядковые номера белковых последовательностей в очереди использования скриптом;
Пункты 3 и 4 повторяются для всех 3-х белков.

Пример запуска скрипта:

$ python3 script2.py X5.fasta P13382 P53942 Q04437

Файл со сборкой X5.fasta я взял отсюда: /P/y19/term3/block2/X5.fasta.

Также прилагаю полученные выходы для всех белков:

# TBLASTN 2.2.28+
# Query: DPOA_YEAST P13382 DNA polymerase alpha catalytic subunit A (2.7.7.7) (DNA polymerase I subunit A) (DNA polymerase alpha:primase complex p180 subunit) (DNA polymerase-primase complex p180 subunit) (Pol alpha-primase complex p180 subunit)
# Database: genome_db
# Fields: query id, subject id, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 7 hits found
DPOA_YEAST	scaffold-423	35.38	1204	641	30	347	1464	582548	579090	9e-176	 589
DPOA_YEAST	scaffold-424	34.77	1208	644	30	347	1464	75468	72007	1e-168	 567
DPOA_YEAST	scaffold-81	30.80	435	251	12	825	1229	342086	340842	6e-40	 162
DPOA_YEAST	scaffold-359	30.57	435	252	12	825	1229	104169	105413	1e-39	 161
DPOA_YEAST	unplaced-816	24.67	300	176	10	847	1101	26390	27274	3e-07	55.1
DPOA_YEAST	scaffold-444	33.87	62	35	2	508	564	571712	571894	1.9	32.7
DPOA_YEAST	scaffold-277	31.58	38	26	0	972	1009	275353	275240	3.8	31.6
# BLAST processed 1 queries

# TBLASTN 2.2.28+
# Query: RNH2A_YEAST P53942 Ribonuclease H2 subunit A (RNase H2 subunit A) (3.1.26.4) (RNase H(201)) (RNase H(35)) (Ribonuclease HI large subunit) (RNase HI large subunit) (Ribonuclease HI subunit A)
# Database: genome_db
# Fields: query id, subject id, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 4 hits found
RNH2A_YEAST	scaffold-199	40.56	143	69	6	152	284	769099	768689	1e-33	93.2
RNH2A_YEAST	scaffold-199	36.67	120	65	4	33	152	769490	769164	1e-33	68.2
RNH2A_YEAST	scaffold-100	39.72	141	73	5	152	284	683313	682903	1e-33	95.5
RNH2A_YEAST	scaffold-100	37.84	111	60	3	42	152	683684	683379	1e-33	65.9
# BLAST processed 1 queries

# TBLASTN 2.2.28+
# Query: KU80_YEAST Q04437 ATP-dependent DNA helicase II subunit 2 (3.6.4.12) (ATP-dependent DNA helicase II subunit Ku80) (High affinity DNA-binding factor subunit 2) (Yeast Ku80)
# Database: genome_db
# Fields: query id, subject id, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 2 hits found
KU80_YEAST	scaffold-423	33.33	57	34	1	243	295	1038260	1038090	1.3	31.6
KU80_YEAST	unplaced-746	24.64	138	82	3	302	422	41434	41036	6.1	29.6
# BLAST processed 1 queries

Также их можно найти по ссылкам: output1.tsv, output2.tsv, output3.tsv.

Исходя из полученных результатов можно судить, что гомологи подобранных мной субъединиц ДНК-полимеразы и РНКазы Н в геноме Amoeboaphelidium protococcarum определённо закодированы, так как мы можем наблюдать довольно высокий, приемлемый процент идентичности и низкие E-value у первых находок в обоих сеансах поиска. В это же время находки по хеликазе II не дали никаких плодотворных результатов: хоть и проценты покрытия не такие низкие, чтобы сходу говорить об отсутствии гомологии, высокий E-value даёт нам окончательно понять, что гомологов этой субъединице хеликазы II в геноме Amoeboaphelidium protococcarum нет.