Упражнения по EMBOSS

Выполненные задания лежат в директории ~/term3/block2/pr9_emboss.txt (адрес относительно домашней папки). Но для удобства продублировал в public_html: Ссылка

Скрипт для решения задачи

Я выбрал третью задачу:
Найти частоты динуклеотидов в геноме бактерии, сравнить их с ожидаемыми и определить динуклеотид, частота которого наиболее отклоняется от наблюдаемой. Подсказка: ожидаемая частота XY = (наблюдаемая частота X) * (наблюдаемая частота Y)
Ссылка на скрипт
Пример выполнения: 
 dimabosov@kodomo:~/public_html/term3$ ./script.bash test.fasta
Count and extract unique words in molecular sequence(s)
Count and extract unique words in molecular sequence(s)
31220
1647 2762 1294 -2063 -2262 -23 598 -1388 -1475 -604 3630 -595 21 -676 -90 -771
-2262 -2063 -1475 -1388 -771 -676 -604 -595 -90 -23 21 598 1294 1647 2762 3630
CG
Описание его работы:
Изначально необходимо ввести имя файла, для которого необходимо выполнить задачу
После для этого файла выполняются 2 команды из EMBOSS, которые записывают количество
динуклеотидов и нуклеотидов в 2 файла: 
    wordcount $filename testt.txt -wordsize=2
    wordcount $filename testtt.txt -wordsize=1
Для каждого нуклеотида была выполнена следующая команда: 
    A=$(($(grep 'A' testtt.txt|awk '{print$2}')+(0)))
$A равняется числу нуклеотидов каждого типа. Если в файле не было данного нуклеотида, то $A не приравнивалось к нулю, а было равно пустому множеству и с ним нельзя было проводить арифметические операции. Именно поэтому я добавлял ноль.
Потом извлекалось число динуклеотидов.
    AT=$(($(grep 'TA' testt.txt|awk '{print$2}')+0))
И теперь считаю их "осклонение" от среднего:
    AT1=$((((($A)*($T)*100000/$((($summ3)))/($summ2)))-((($AT)*100000/($dinucl)))))
Домножаю на 100000 из-за того, что bash не умеет работать с дробными числами. Я попытался найти способ работы с дробными числами и он достаточно запутанный. Посчитал этот метод гораздо удобнее, ибо если я увеличу отклонение от ожидаемого у каждого динуклеотида на одно и тоже число, то самое большое значение всё ещё будет самым большим, и я его не потеряю.
При подсчёте частоты динуклетидов у второго нуклеотида частота выражается как: "число нуклеотидов/общее число нуклетиодов", а общее число нуклеотидов на 1 меньше, так как мы уже знаем первый. При больших последовательностях это почти не ощущается, но если запустить, например, последовательность из 10 букв, то становится заметно.
    otcl=($AT1 $TA1 $AC1 $CA1 $AG1 $GA1 $GT1 $TG1 $CT1 $TC1 $CG1 $GC1 $GG1 $AA1 $CC1 $TT1)
    otcl1=(AT TA AC CA AG GA GT TG CT TC CG GC GG AA CC TT)
Вводим 2 массива: значения отклонений и динуклеотиды, так чтобы они соотносились между собой(k-тому отклонению первого массива соответствовал k-тый динуклеотид.) 
    IFS=$'\n' sorted=($(sort -n <<<"${otcl[*]}")); unset IFS
Сортируем список из всех отклонений.
И теперь находим наибольшое и наименьшее отклонение. Домножаем наименьшее на -1 и сравниваем их, при этом определив их положние в неотсортированном списке с отклонениями. Мы ищем его положение в неотсортированном списке из-за того, что при сортировке теряется соответсвие динуклеотид- его частота. И найдя самое большое значение, выводим динуклотид их списка.
    rm testt.txt
    rm testtt.txt
Удалим 2 файла, полученные в процессе. При тестированнии программы заметил, что если ввести несуществующий файл, то скрипт выдаст, что он его не нашёл, но при этом выдаст вполне полноценный результат(динуклетид). Это было из-за того, что раньше эти 2 файла сохранялись и ей было с чем работать, даже если его не перезаписывали. Теперь же она выдаст ошибку и даст явно понять, что введен несуществующий файл. Можно было бы подсчёт всех частот и сумму завести под цикл, но при применении скрипта на последовательности длиной 50-100 тысяч основании, результат выдается моментально, поэтому посчитал, что можно оставить так.

BLAST+

Необходимо было найти гомологи белка в неаннотированной сборке генома Amoeboaphelidium protococcarum. Так как это гриб (вернее родственник) было решено провести сравнение с белками хорошо изученного Saccharomyces cerevisiae. Скачаем несколько белков из uniprot c помощью команды:
    seqret sw:AC_PROTEIN NAME_PROTEIN.fasta
Аконитаза
Актин
Гексокиназа-2
Создадим локальную базу данных: 
    dimabosov@kodomo:~/public_html/term3/pr9$ makeblastdb -in X5.fasta -dbtype nucl

    Building a new DB, current time: 11/09/2020 15:53:45
    New DB name:   X5.fasta
    New DB title:  X5.fasta
    Sequence type: Nucleotide
    Keep Linkouts: T
    Keep MBits: T
    Maximum file size: 1000000000B
    BLAST options error: File X5.fasta does not exist
Ну и теперь выполним сравнение белков с геномом, использовав при этом tblastn: 
   tblastn -query NAME_PROTEIN.fasta -db X5.fasta > NAME_PROTEIN.txt
Файлы, полученные этим методом:
Аконитаза
Актин
Гексокиназа-2
Чтобы сделать вывод о том, что эти белки являются гомологами, достаточно можно посмотреть на их E-value.
0.0 , 0.0 и 1e-96 для трёх белков соответственно.Помимо этого нужно ещё посмотреть на процент идентичных аминокислот: 68, 89 и 40 соответственно. Достаточно сложно сказать о готологичности у гексокиназы-2, но всё же, столь большой e-value и то, что похожих по свойствам аминокислот при выравнивании (параметр Positives) равен 60% может говорить о том, что этот белок у Amoeboaphelidium protococcarum претерпел большое количество мутаций и есть шанс, что он является гомологом. Про остальные 2 белка. Думается, что они могут иметь гомологи в геноме Amoeboaphelidium protococcarum с большой вероятностью. К тому же у всех 3 белков выранивание происходило почти по всей длине,а это очень хорошо.