Выбранные упражнения по применению EMBOSS доступны в файле ~/term3/block2/pr9_emboss.txt, а также здесь.
Выбрана третья задача: найти частоты динуклеотидов в геноме бактерии, сравнить их с ожидаемыми и определить динуклеотид, частота которого наиболее отклоняется от наблюдаемой.
Скрипт доступен по ссылке
Команда запуска:
python3 dinucleotide_freq.py [sequence]
Последовательность задаётся в формате USA с путём к файлу (если чтение из файла) относительно расположения скрипта. Выдача выглядит так:
Most derivating dinucleotide: CG
Если добавить параметр -verbose (на самом деле, если добавить после последовательности что угодно), то помимо динуклеотида с наибольшим отклонением выведутся также таблицы частот, ожидаемых частот и отклонения от ожидаемого для всех динуклеотидов:
Frequences:
[[1 \ 2 A T G C ]
[A 0.0794 0.0604 0.0747 0.0509 ]
[T 0.0472 0.0792 0.0768 0.0612 ]
[G 0.0626 0.0509 0.0681 0.0542 ]
[C 0.0762 0.0738 0.0163 0.0668 ]]
Exspected:
[[1 \ 2 A T G C ]
[A 0.0704 0.0702 0.0626 0.0619 ]
[T 0.0702 0.0699 0.0624 0.0616 ]
[G 0.0626 0.0624 0.0556 0.055 ]
[C 0.0619 0.0616 0.055 0.0543 ]]
Difference:
[[1 \ 2 A T G C ]
[A 0.0089 0.0097 0.0121 0.0109 ]
[T 0.023 0.0093 0.0144 0.0005 ]
[G 0.0 0.0114 0.0125 0.0008 ]
[C 0.0143 0.0122 0.0387 0.0125 ]]
Most derivating dinucleotide: CG
Пример выше - выдача сборки 22 хромосомы человека (/P/y19/term3/block2/pr9/hs_ref_GRCh38.p7_chr22.fa).
Применение программы к простой тестовой последовательности ATTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT даёт вполне согласующийся с логикой результат:
Frequences:
[[1 \ 2 A T G C ]
[A 0.0 0.0238 0.0 0.0 ]
[T 0.0 0.9524 0.0 0.0 ]
[G 0.0 0.0 0.0 0.0 ]
[C 0.0 0.0 0.0 0.0 ]]
Exspected:
[[1 \ 2 A T G C ]
[A 0.0006 0.0232 0.0 0.0 ]
[T 0.0232 0.9529 0.0 0.0 ]
[G 0.0 0.0 0.0 0.0 ]
[C 0.0 0.0 0.0 0.0 ]]
Difference:
[[1 \ 2 A T G C ]
[A 0.0006 0.0006 0.0 0.0 ]
[T 0.0232 0.0006 0.0 0.0 ]
[G 0.0 0.0 0.0 0.0 ]
[C 0.0 0.0 0.0 0.0 ]]
Most derivating dinucleotide: TA
Действительно, нуклеотиды G и C последовательности отсутствуют, соответственно, их ожидаемые и наблюдаемые частоты равны нулю. Из динуклеотидов, составленных из А и Т отсутствуют АА и ТА, однако, поскольку Т имеет большую частотность, чем А, ожидаемая частотность ТА больше, чем АА, соотвественно, разница с нулём больше и именно он опознаётся как наиболее отклоняющийся.
Из сборки генома Amoeboaphelidium protococcarum была сформирована база, по которой далее производился BLAST со следующими параметрами:
tblastn -query RPB1_YEAST.fasta -db X5.fasta -out RPB1_YEAST.out -evalue 1 -outfmt '7 qseqid sseqid evalue qstart qend qcovs sstart send slen' -num_alignments 1000000
Поиск белков производился по Uniprot c повыщением уровня таксономии, если на текущем не было белков из swissprot. В итоге все белки найдены для Opisthokonta. Далее последовательность получалась командой entret "fasta::swissprot:PA1_HUMAN" PA1_HUMAN.fasta.
Выбраны три белка:
Днк-полимераза альфа (субъединица А) из Saccharomyces cerevisiae. Последовательность, выдача. Присутствует выравнивание с хорошим (9e-176) e-value и достаточно большим (76%) покрытием, скорее всего, это гомолог (два, на 423 и 424 скэффолдах). Кусочки с покрытием 28%, возможно, имеют просто гомологичный домен. Видно, что выравнивания не находятся на краю скэффолда и не продолжают друг друга.
Днк-зависимая РНК полимераза II (субъединица RPB1) из Saccharomyces cerevisiae. Последовательность, выдача. Присутствует выравнивание с хорошим (0.0) e-value и достаточно большим (88%) покрытием, скорее всего, это гомолог (два со схожими значениями, на 300 и 157 скэффолде, остальные хуже).
Днк-зависимая РНК полимераза III (субъединица RPC) из Homo sapiens. Последовательность, выдача. Присутствует выравнивание с хорошим (0.0) e-value и 99% покрытием. Снова два гомолога с одинаковыми значениями (остальных хуже), на тех же 300 и 157 скэффолдах (но, заметим, что координаты выравниваний по subject не пересекаются, то есть, это не гомология двух РНК-полимераз с одним участком хромосомы).