Для проведения опыта проводилось выравнивание маленького протеома по большому. В качестве маленького изначально использовался полный протеом Haemophilus influenzae, но из-за трудоемкости задачи, были задействованы только первые 10% строк.
Для выделения этих строк использовались стандартные функции bash - wc и head.
wc -l haemofilusproteom.fasta 11336 haemofilusproteom.fasta
head -1133 haemofilusproteom.fasta > h10.fasta
Таким образом были получены первые 10% (1133) строк. Далее с помощью текстового редактора vi, последняя fasta последовательность была вырезана, чтобы все записи были целыми.
Для индексирования был использован протеом крысы (Rattus norvegicus). Индексирование проводилось с помощью команды makeblastdb.
makeblastdb -in ratproteom.fasta -dbtype prot -out indrat
Все выравнивания отбирали результаты с E-value не более 0.001 (опция -evalue 0.001) и записывали результат в табличном формате (опция -outfmt 6), при этом фиксируя время выполнения с помощью 'time'. Для подсчета количества последовательностей использовалась команда 'wc -l'.
Первое выравнивание было эталонным, никакие другие параметры при его подсчете изменены не были. Команда для расчета представлена непосредственно ниже, а табличный файл можно просмотреть по ссылке, в нем 1098 строки. На построение этих выравнивании затрачено около 45 с (полный файл с временами исполнения находится по адресу).
time blastp -query h10.fasta -db indrat -evalue 0.001 -outfmt 6 -out t0
Далее представлены выавнивания с разным размером слова (опция -word_size int). Всего было получено 4 серии выравниваний: со значением параметра 7,6,3 (что соовтетствует значению по-умолчанию) и 2. Количество строк в этих файлах равняется 1105, 1105, 1098, 1105 соответственно, а времена работы - 9 мин 28 с, 11 мин 36 с, 45 с, 1 мин 47 с. Типичная команда выглядит как:
time blastp -query h10.fasta -db indrat -evalue 0.001 -word_size 7 -outfmt 6 -out t1
Следующим параметром стал мининмальный порог сходства слова threshold (опция -threshold int изначально 11). Всего были измерены серии вырваниваний со значениями threshold=10,1 и 15. В них были найдены 1105, 1105 и 885 выравниваний соответственно, на это было затрачено 1 мин 31с, 161 мин 41 с и 8 с. Типичная команда выглядит как:
time blastp -query h10.fasta -db indrat -evalue 0.001 -thresold 10 -outfmt 6 -out t5
Также был рассмотрен параметр максимального расстояния между сходными словами window size (опция -window_size int изначально 40). Были замерены количества выравниваний для знаений window size=50,10,20,40,60 и 150, им соответствуют времена 48 с, 14 с, 36 с, 45 с , 51 с, 1 мин 14 с и количества строк 1099, 1055, 1090, 1098, 1099, 1099. Первая из ряда команд:
time blastp -query h10.fasta -db indrat -evalue 0.001 -window_size 50 -outfmt 6 -out t8
Следующий параметр - порог отклонения от максимального веса (в битах) при расширении от якоря xdrop ungap (опция -xdrop_ungap float изначально 7). Рассмотренные значения параметра: 1,5,10,15 и 25. Затраченное время - 31 с, 42 с, 47 с, 50 с , 53 с. Количества найденных выравниваний - 933, 1096, 1098, 1098, 1098. Типовая команда:
time blastp -query h10.fasta -db indrat -evalue 0.001 -xdrop_ungap 1 -outfmt 6 -out t14
Далее были рассмотрен параметр предварительного отклонения от веса при расширении с гэпами xdrop gap (опция -xdrop_gap float изначально 15 бит). Использованные значения: 1,5,10,15 и 25, времена исполнения: 37 c, 37 c, 40 c, 45 c и 58, количества строк: 544, 544, 936, 1098 и 1133. Типичная команда выглядит как.
time blastp -query h10.fasta -db indrat -evalue 0.001 -xdrop_gap 1 -outfmt 6 -out t19
Помимо xdrop gap были произведены запуски blastp с различными значениями параметра xdrop gap final (опция -xdrop_gap_final изначально 25). Предположительно, этот параметр устанавливает дополнительный порог отклонения для веса (в битах) уже после того, как сработал предыдущий, а построение выравнивания закончено, и позволяет "отсеять" еще некоторые выравнивания. Предположение основано на описании параметров при запуске функции для PSI-BLAST по адресу. Точно установлено, что количества последовательностей при одинаковых значениях параметров xdrop gap и xdrop gap final различаются. Так, при запусках со значениями второго, равными 1,5,10,15,16.46875 (число являлось предположительным дефолтным значением при условии монотонности функции, что впоследствии оказалось неверным), 20 и 25 количества выравниваний составляют 525, 553, 889, 1071, 1098, 1104, 1098 соответственно. Этот параметр в немалой степени степени повлиял на количество найденных выравниваний уменьшив его более, чем в 2 раза по сравнению с запросом с неизмененными параметрами. Несколько вразрез с предположением идет тот факт, что количество выданных строк при одинаковых значениях 5 и 10 оказалось больше, чем у xdrop gap, а ожидалось меньшее или равное. Помимо этого, при значении 30 количество строк оказалось меньше, чем при 25, хотя ожидалось, что при увеличении xdrop gap final количество найденных выравниваний будет возрастать, так как отклонение будет расти и позволит выбрать новые варианты. Изменение параметра почти не повлияло на время работы, - оно изменялось в пределах 1 с и составило 44 с, 44 с, 44 с, 45 с, 45 с, 45 с, 45 с. При этом между 15 и 16.46875 даже был небольшой спад в 43 мс, вероятно, погрешность. Одна из выполненных команд:
time blastp -query h10.fasta -db indrat -evalue 0.001 -xdrop_gap_final 1 -outfmt 6 -out t24
Далее рассмотатривался параметр, учитывающий сложность состава белков (и low complexity regions; опция -comp_based_stats compo по-умолчанию conditional compositional score matrix adjustment = 2). В Blast от NCBI этот параметр можно установить равным одному из четырех вариантов: no adjustment, composition based statistics, conditional compositional score matrix adjustment или universal compositional score matrix adjustment, все они были использованы. При этом они все дали одинаковый результат - 1676 строк (рекордно высокое) и округленные до целых 44 с работы, в том числе запуск с третьем значением параметра. Возможно, compo - не просто строка, и отсутствие ошибки при выполнении не говорит о правильном введении параметра. Все четыре команды представлены ниже:
time blastp -query h10.fasta -db indrat -evalue 0.001 -comp_based_stats 'universal compositional score matrix adjustment' -outfmt 6 -out t34
time blastp -query h10.fasta -db indrat -evalue 0.001 -comp_based_stats 'conditional compositional score matrix adjustment' -outfmt 6 -out t33
time blastp -query h10.fasta -db indrat -evalue 0.001 -comp_based_stats 'composition-based statistics' -outfmt 6 -out t32
time blastp -query h10.fasta -db indrat -evalue 0.001 -comp_based_stats 'no adjustment' -outfmt 6 -out t31
Кроме того, были использованы разные матрицы начисления весов (опция -matrix matrix_name изначально BLOSUM62): PAB30, PAM70, PAM250, BLOSUM80, BLOSUM62, BLOSUM45, BLOSUM50 И BLOSUM90. Они показали 294, 556, 728, 919, 1098, 1014, 1084 и 884 выравнивания соответственно и потратили 9 с, 17 с, 2 мин 18 с, 26 с, 45 с, 36 с, 4 мин 7 с и 1 мин 20 с. Эти результаты эмпирически подтвердили выбор BLOSUM62 в качестве лучшей матрицы - она обнаружила больше всех выравниваний и при этом потратила удовлетворительное время. Если требуется быстрота выполнения можно воспользоваться PAM30. Типичная команда:
time blastp -query h10.fasta -db indrat -evalue 0.001 -matrix 'PAM30' -outfmt 6 -out t35
Предпоследний рассматриваемый параметр - штраф за открытие гэпа (опция -gapopen open_penalty изначально 11). Интересно что он и штраф за расширение гэпа не могут принимать любые значения, за ними закреплены конкретные пары для каждой матрицы. Например, для BLOSUM62 пара (1;1) не является допустимой:
BLAST query/options error: Gap existence and extension values of 1 and 1 not supported for BLOSUM62 supported values are: 32767, 32767 11, 2 10, 2 9, 2 8, 2 7, 2 6, 2 13, 1 12, 1 11, 1 10, 1 9, 1
Ввиду ограниченности значений параметра были рассмотрены 9, 10, 11, 12 и 13, им соответствуют количества строк 1101, 1131, 1098, 1063 и 1039, а также времена исполнения 52 с, 48 с, 46 с, 44 с и 43 с. Вполне ожидаемо, что с ростом штрафа уменьшается количество выравниваний. Типичная команда выглядит как:
time blastp -query h10.fasta -db indrat -evalue 0.001 -gapopen 9 -outfmt 6 -out t43
Последним параметром стал штраф за расширение гэпов (опция -gapextend extend_penalty изначально 1), он связан с предыдущим как идейно (расширение невозможно без начала), так и через фиксированные значения. Все были рассмотрены gap extend penalty=1, и 32767 одновременно с предыдущим параметром, были получены 1098, 1024 и 997 строк за 46 с, 41 с и 38 с. Последнее значение затронуло два измененных параметра и установило настолько крупный штраф, что, вероятно, все выравнивания построены без гэпов. Типичная команда для изменения параметра:
time blastp -query h10.fasta -db indrat -evalue 0.001 -gapextend 2 -outfmt 6 -out t49
На самом деле -comp_based_stats принимает значения 0, 1, 2, 3, слегка поздно был обнаружен мануал по параметрам BLAST+ от NCBI. Поэтому были проведены дополнительные расчеты, все четыре значения параметра показали одинаковое время исполнения 45 с, количетсва строк в выдачах равны соответственно 1676, 1139, 1098 и 1100.
Кроме того, для blastn значение параметра xdrop gap final по-умолчанию равняется 100, что подтолкнуло меня к рассмотрению значений 50,100 и 200. С увеличением параметра немного растет время работы: 47 с, 49 с, 53 с, количество строк почти не меняется: 1089, 1084, 1084. Описание в указанном мануале совпадает с приведенным выше - "эвристическое значение (в битах) для окончательного выравнивания с гэпами", но почему-то для blastp не указан параметр xdrop gap, только xdrop gap final. Возможно, эти параметры вместе имеют смысл лишь при выравниваний нуклеиновых кислот (blastn).
На время работы влияют не так много параметров: word size, threshold (со значением 1 работает долго, потому что рассмтаривает много вариантов, даже с 10 % протеома заняло больше 2.5 часов) и matrix (что неожиданно, потому что сложность алгоритма не должна значительно изменяться). При этом похоже, что BLAST не рассматривает некоторые выравнивания при уменьшении word size, так как при уменьшении этого параметра время работы должно было монотонно увеличиваться, но этого не происходило. На число найденных выравниваний значительно повлияли параметры: xdrop gap, xdrop gap final и matrix (чтобы сильно понизить число еще можно использовать word size=1, но частный случай, вероятно остальные параметры тоже так могут). Точное различие между xdrop gap и xdrop gap final установлено не было, но оно гарантировано существует.
В качестве банка данных использовался уже проиндексированный протеом крысы, а в качестве изначальной последовательности - первичная структура уридин фосфорилазы из холерного вибриона (Vibrio cholerae), её длина 253 аминокислотных остатка, что близко к среднему значению среди всех белков. Для подготовки серии опытов был скачана последовательность белка:
seqret uniprot:q9k4u1 uf -auto
Далее последовательность была случайным образом перемешана 100 раз и записана в другой файл.
shuffleseq uf pp -shuffle 1000
В связи с рекомендациями не использовать функцию 'seqretsplit', была использована команда 'seqret' с опцией -ossingle 1. Но поскольку для второй команды fasta-последовательности с одним названием считаются одинаковыми, потребовалось сделать имена уникальными, например, пронумеровать с использованием Python. В результате образовался файл pp1. Уже к нему была применена выше упомянутая команда:
seqret pp1 a -ossingle 1
Это образовало 1000 файлов с типовым названием q9k4u1_vibcli.fasta, где i - номер сгенерированной последовательности.
Далее 1000 раз был запущен BLAST с помощью Bash-скрипта (представлена лишь шестая версия скрипта с порогом E-value 8; четвертая строка в нем - просто отладка; в процессе написания также были записаны 1000 файлов с выравниваниями для порога 0.25, но так как последовательности бессмыслены я не стал их указывать, они носят имена ti, где i - число от 50 до 1050, соответсвующее i-50 сгенерированной последовательности). Каждый раз скрипт кроме отладки выдает количество находок, оно составило 45, 109, 245, 480, 1012 и 1892 соответственно порогам E-value 0.25, 0.5, 1, 2, 4 и 8. Зная, что в каждом опыте принимало участие 1000 последовательностей, среднее число находок можно найти как частное суммы всех находок за проход скрипта и 1000. В таком случае средние количества находок равны: 0.045, 0.109, 0.245, 0.48, 1.012 и 1.892.
При сравнении этих значений с порогами E-value можно отметить, что пороги примерно в 4-5 раз больше среднего количества находок в каждом случае. Это можно объяснить отличием проиндексированного крысинного протеома от идеального банка, в котором бесконечное количество независимых случайных последовательностей. Примечательно, что отличие в 4-5 раз по сравнению с экспоненциальными значениями E-value влияет несущественно (всего лишь O(10k)), а это значит, что протеом крысы - относительно неплохой банк.