Феодальный кризис сыграл на руку слонам,
Превратив Афины в безнадёжный Суринам.
Кобыла и трупоглазые жабы, "Гудзон (Вспоминая Варанов)"

Написание скрипта, который что-то делает

 Была выбрана следующая задача: «Проверить, сколько находок с E-value < 0.1 в среднем находит blastn для случайной последовательности данной длины в данном геноме бактерии».
 Я решил писать скрипт на bash, потому подумал, что научиться считать числа в bash будет проще, чем научиться вызывать другие программы в Python, но я ошибся. Тем не менее, скрипт я написал. Посмотреть на скрипт в браузере можно по этой ссылке, а скачать – по этой.
Скрипт работает следующим образом:

1. Генерируется 100 случайных последовательностей заданной длины.
2. Каждая последовательностей бластится по заданному геному бактерии.
3. Количество подходящих находок складывается и делится на 100.

Первым аргументом скрипт принимает длину последовательности, а второй файл с геномной последовательностью. Пример запуска:

legoushque@kodomo:~/public_html/term3/pr9/1$ bash script.bash 100 Rickettsia_rickettsii.fasta
The average number of alignmenst with E-value less than 0.1: .08

legoushque@kodomo:~/public_html/term3/pr9/1$ bash script.bash 10000 Rickettsia_rickettsii.fasta
The average number of alignmenst with E-value less than 0.1: .12 

Видно, что в обоих случаях число находок меньше единицы, а значит в большинстве случаев находок с E-value меньше 0.1 вообще нет.

Поиск гомологов белков в неаннотированном геноме

 Стояла задача найти в геноме Amoeboaphelidium protococcarum гомологов каким-нибудь трём белкам, которые должны быть у всех эукариот. Для начала нужно было установить, из каких организмов вообще брать белки. Поскольку NCBI Taxonomy не даёт никаких таксонов между Aphelida и Opisthokonta было принято решение провести википедийный анализ ситуации (Рис. 1). Было решено брать белки из Eumycota.

 Были выбраны следующие белки (в скобках поисковый запрос в UniProt, по которому нужный белок был найден): актин("actin fungi"), бета цепь тубулина("tubulin fungi") и бета субъединица АТФ-синтазы ("atp synthase fungi"). Так получилось, что все белки из Saccharomyces cerevisiae. Последовательности скачивались следующим образом:

legoushque@kodomo:~/public_html/term3/pr9/2$ seqret sw:P60010 proteins/P60010.fasta
Read and write (return) sequences
legoushque@kodomo:~/public_html/term3/pr9/2$ seqret sw:P02557 proteins/P02557.fasta
Read and write (return) sequences
legoushque@kodomo:~/public_html/term3/pr9/2$ seqret sw:P00830 proteins/P00830.fasta
Read and write (return) sequences

 Ссылки на последовательности:

1. Актин
2. Бета цепь тубулина
3. Бета субъединица АТФ-синтазы

 Далее индексируем сборку генома, в которой будем искать гомологи:

legoushque@kodomo:~/public_html/term3/pr9/2$ makeblastdb -in X5.fasta -dbtype nucl


Building a new DB, current time: 11/08/2020 22:15:12
New DB name:   X5.fasta
New DB title:  X5.fasta
Sequence type: Nucleotide
Keep Linkouts: T
Keep MBits: T
Maximum file size: 1000000000B
Adding sequences from FASTA; added 1868 sequences in 0.85028 seconds.

 Поскольку мы ищем по белковой последовательности в нуклеотидной базе данных, будем использовать tblastn:

legoushque@kodomo:~/public_html/term3/pr9/2$ for x in `ls proteins`; do tblastn -query $x -db X5.fasta > $x.txt; done  

 Выдачи можно увидеть по ссылкам ниже:

1. Актин
2. Бета цепь тубулина
3. Бета субъединица АТФ-синтазы

 Для каждого из белков есть находки с E-value < e-100 практически на всю длину, что говорит о том, что для каждого из этих белков в геноме Amoeboaphelidium protococcarum есть гомологи.