Работа в EMBOSS
Упражнения по EMBOSS
Необходимо было сделать 7 из 16 упражнений, представленных в задании. Файл с ними недоступен для чтения людям без доступа к серверу, поэтому я прилагаю копию ниже:
2 secretsplit input.fasta -auto 4 transeq input.fasta output.fasta -table5 getorf input.fasta output.fasta -minimize 8 featcopy gb::input.gb ggf::output.gff 9 extractfeat input.gb output.fasta -type exon 11 makenucseq -amount 3 -length 100 -outseq output.fasta -auto 12 cusp input.fasta output.cusp
Скрипт в bash, решающий определённую задачу
Я выбрал первое задание со следующей формулировкой:
1. Проверить, сколько находок с E-value < 0.1 в среднем находит blastn для случайной последовательности данной длины в данном геноме бактерии. (Подсказка: сгенерируйте 100 случайных последовательноcтей и задайте их как query, выдачу blastn задайте табличную).
В результате на языке bash был написан скрипт evrgblastn.bash, принимающий первым аргументом количество последовательностей, которое должно быть в созданном случайным образом банке, а вторым длины случайных последовательностей в банке. Третьим аргументом задаётся запись с последовательностями из генома в формате USA (в скрипте используются команды EMBOSS). На выходе последним сообщением выводится среднее количество (до двух знаков после запятой) находок в банке из заданного с E-value не больше 0.1 для случайной последовательности с заданной длиной. Стоит отметить, что в ходе работы создаются новые файлы (в том числе база данных из генома), автоматически удаляющиеся в конце работы. Ниже приведён пример использования скрипта с заданным аргументами два раза подряд (каждый раз будет создано 100 последовательностей с длиной 100 нуклеотидов):
1 запуск:
starlingsden@kodomo:~/term3/block3$ ./evrgblastn.bash 100 100 ddbj:CP027599
Building a new DB, current time: 12/22/2020 00:56:11
New DB name: outseq.fasta
New DB title: outseq.fasta
Sequence type: Nucleotide
Keep Linkouts: T
Keep MBits: T
Maximum file size: 1000000000B
Adding sequences from FASTA; added 1 sequences in 0.108709 seconds.
Среднее количество находок с E-value не больше 0.1 при данных параметрах равно 0
2 запуск:
starlingsden@kodomo:~/term3/block3$ ./evrgblastn.bash 100 100 ddbj:CP027599
Building a new DB, current time: 12/22/2020 00:45:06
New DB name: outseq.fasta
New DB title: outseq.fasta
Sequence type: Nucleotide
Keep Linkouts: T
Keep MBits: T
Maximum file size: 1000000000B
Adding sequences from FASTA; added 1 sequences in 0.108346 seconds.
Среднее количество находок с E-value не больше 0.1 при данных параметрах равно .02
Необходимо учитывать, что точность измерений растёт при увеличении количества последовательностей. Так, за несколько осуществлённых мной запусков для 500 последовательностей ни разу не было значения, отличающегося от 0, в то время как для 100 последовательностей довольно часто возникали значения от 0.02 до 0.09.
Скрипт можно найти в директории /home/students/y19/tarlingsden/term3/block3, но к нему нет доступа у людей без возможности зайти на сервер kodomo, поэтому я дополнительно прилагаю ссылку на файл с кодом: evrgblastn.bash.
Работа с BLAST+
Поиск гомологов белков в неаннотированном геноме
Необходимо было провести поиск гомологов белков в неаннотированной сборке генома Amoeboaphelidium protococcarum.
Так как данный организм является родственником грибов, было решено искать достаточно хорошо аннотированных организмов для поиска гомологов их белков среди царства Fungi (таксон может не соответствовать современным представлениям о систематике). Для работы был выбран Encephalitozoon cuniculi.
После запроса "taxonomy:"Encephalitozoon cuniculi (Microsporidian parasite) [6035]" AND reviewed:yes" на сайте UniProt появилась таблица с аннотированными белками. Среди них достаточно важными, чтобы почти наверняка присутствовать в геноме исходного организма, показались убиквитин (Q8SWD4) и 40S рибосомальный белок S5 (Q8SS72). Также к ним добавляется каталитическая субъединица А протонной АТФ-азы V-типа (Q8SQU9).
Последовательности соответствующих белков были скачаны с ипользованием команды (для последнего белка) entret "uniprot:Q8SQU9" stdout 2> /dev/null | seqret stdin Q8SQU9.fasta 2> /dev/null (2> /dev/null удаляет сообщения EMBOSS, которые выводятся на StdErr).
Для использования BLAST+ необходимо было сначала сделать из сборки генома бактерии X5.fasta базу данных с помощью команды makeblastdb -in X5.fasta -dbtype nucl, после чего можно было преступить к бласту. Т.к. мы сравниваем последовательности исходно белковые с исходно нуклеотидными, использовалась команда tblastn: tblastn -query Q8SWD4.fasta -db X5.fasta -db_gencode 6 > Q8SWD4.txt (пример для убиквитина, "-db_gencode 6" задаёт соответствующую Amoeboaphelidium protococcarum таблицу генетического кода).
1. Убиквитин
В выдаче tblastn представлены примерно одинаковы епо характеристикам выравнивания на всю длину убиквитина (1-76), где невооружённым глазом видно почти полное соответствие. Вот один из примеров:
> scaffold-105
Length=655906
Score = 150 bits (378), Expect = 1e-42, Method: Compositional matrix adjust.
Identities = 74/76 (97%), Positives = 75/76 (99%), Gaps = 0/76 (0%)
Frame = -3
Query 1 MQIFVKTLTGKTITLEVEPSDSIENVKAKIQDKEGIPPDQQRLIFAGKQLEDGRTLSDYN 60
MQIFVKTLTGKTITLEVE SD+IENVKAKIQDKEGIPPDQQRLIFAGKQLEDGRTLSDYN
Sbjct 655664 MQIFVKTLTGKTITLEVESSDTIENVKAKIQDKEGIPPDQQRLIFAGKQLEDGRTLSDYN 655485
Query 61 IQKESTLHLVLRLRGG 76
IQKESTLHLVLRLRGG
Sbjct 655484 IQKESTLHLVLRLRGG 655437
Гомология очевидна, почти полная идентичесность позволяет утверждать, что найденная в сборке генома последовательность соответствует убиквитину.
2. 40S рибосомальный белок S5
Для этого белка необходимо учитывать, что уровень аннотации говорит лишь о существовании белка, в последовательности могут содержаться ошибки, что в теории может привести к снижению количества идентичных и схожих аминокислотных остатков.
Обратившись к выдаче tblastn, снова обнаружим выравнивание по всей длине белка (1-189), почти без гэпов (только два), что позволяет предположить сохранение общей структуры. Небольшое значение E-value (4e-45, 5e-45) вместе с достаточно высокой идентичностью с учётом недостаточной аннотированности говорит о том, что гомология скорее всего есть, но, чтобы понять, что найденная в геноме последовательность соответствует 40S рибосомальному белоку S5, необходимо знать о характере его консервативности.
3. Субъединица А протонной АТФ-азы V-типа
Как видно в выдаче tblastn, выравнивание прошло почти по всей длине (22-606 из 607), при этом для такой большой последовательности идентичность составляет 65% при машинном нуле у E-value. Также можно обратить внимание на крайне низкую частоту гэпов (всего 2) и почти полное соответствие последовательностей друг другу ближе к середине. Из всего это следует, что перед нами гомологи. Также почти наверняка мы нашлии таким образом в геноме последовательность, кодирующую субъединицу А протонной АТФ-азы, однако это также необходимо перепроверить, изучив разнообразие АТФ-аз.
