Обновлено 21.08.2012
В статье пойдет речь о реализации простого алгоритма поиска регуляторных мотивов на языке программирования Haskell. Попутно с решением задачи постараюсь также разьяснить некоторые особенности языка.
Регуляторные мотивы представляют собой короткие последовательности нуклеотидов, которе, как правило, расположены непосредственно перед генами и используются для контроля экспрессии генов. Т.е. регуляторные мотивы задают условия, при которых ген будет "включаться" или "выключаться".
Подробно ознакомиться с задачей можно в замечательной книге [2], а также в презентации. Желающих подробнее ознакомиться с языком Haskell отсылаю к книгам [1] и [3].
Задача формулируется следующим образом:
в заданном наборе из m последовательностей ДНК длины n каждая, необходимо найти неизвестную последовательность длиной k (k-мер), которая встречается во всех m последовательностях. При этом k-мер может быть представлен неточно, т.е. допускается наличие мутаций.
Будем решать эту задачу сформулированную в форме поиска median string. Понятие median string тесно связано с расстоянием Хэмминга: это число позиций, в которых соответствующие символы двух слов одинаковой длины различны. Например, расстояние Хэмминга dH(“AAAA”, “AACA” ) = 1. Рассматриваются только строки одинаковой длины.
Таким образом расстояние Хэмминга является мерой различия строк. Теперь определим общее расстояние Хэмминга (total Hamming distance). Для этого будем действовать следющим образом:
Тогда задачу поиска median string можно сформулировать как
найти k-мер, минимизирующий общее расстояние Хэмминга для m последовательностей ДНК длиной n.
В данной статье мы рассмотрим алгоритм поиска мотива методом полного перебора (или “грубой силы”).
Исходя из формулировки задачи мы уже можем наметить вариант решения:
Теперь у нас есть все исходные данные для решения задачи, приступим к написанию программы.
Я буду придерживаться стиля написания программы bottom-up: сначала реализуются функции нижнего уровня, а затем из них, как из строительных блоков собираются функции верхнего уровня.
В первую очередь необходимо посчитать расстояние Хэмминга. Как следует из определения, нам необходимо найти количество позиций в которых две строки различны.
dH seq1 seq2 = sum $ zipWith diff seq1 seq2 where diff a b | a==b = 0 | otherwise = 1
Для построения функции dH используем стандартные функции языка Haskell для работы со списками: zipWith - применяет функцию diff к элементам последовательностей seq1 и seq2, находящимся в одинаковых позициях(почти как zip, только применяет функцию к сформированым парам). Таким образом, формируется список содержащий 0, если элементы равны и 1, если не равны. Далее суммируем элементы полученного списка, что и даст расстояние Хэмминга. Чтобы передать результат zipWith на вход функции sum используется опрератор “$”. Он нужен здесь только для того, чтобы писать меньше скобок, т.е. sum $ zipWith diff seq1 seq2 эквивалентно sum (zipWith diff seq1 seq2)
Чтобы посмотреть, как это выглядит воспользуемся REPL для языка Haskell и посмотрим результат применения функции dH к различным входным данным.
*BruteforceMotifSearch> dH "AAAA" "AAAA" 0 *BruteforceMotifSearch> dH "AAAA" "AAAC" 1 *BruteforceMotifSearch> dH "AAAA" "TAAC" 2
Далее нам нужна функция, умеющая из заданной строки получать все подстроки длиной k:
k_mers k seq = map (take k) $ take (n-k+1) $ tails seq where n = length seq
Разберем последовательно, что здесь происходит. Сначала вызывается функция tails, возвращающая список “хвостов” переданного на вход списка:
*BruteforceMotifSearch> tails "ACGT" ["ACGT","CGT","GT","T",""]
Все “хвосты” не нужны, нас интересуют только только те, которые имеют длину не меньше k. Таких “хвостов” будет n-k+1, где n - длина списка с “хвостами”, поэтому с помощью функции take берем только первые n-k+1 элементов. Но это еще не то, что нам нужно. Теперь воспользуемся функцией map - одной из базовых функций в функциональных языках программирования. Её задача, в общем, довольно проста: получить новый список из исходного, применяя к каждому элементу некоторую функцию, переданную также в качестве параметра. Например:
*BruteforceMotifSearch> map (\x -> x*x) [1,2,3] [1,4,9]
Здесь (\x -> x*x)- анонимная функция, которая считает квадрат переданного ей числа. Применив её к списку чисел получим список квадратов.
В нашем случае с получением всех k-меров нужно взять k первых символов из каждого “хвоста”. Поэтому в качестве функции, которая будет применяться к каждому элементу со списком хвостов будем использовать (take k). Да-да, это именно функция. Это явление в функциональных языках называется каррирование: функция take принимает 2 аргумента - количество элементов и список, - если мы передадим только первый аргумент, то получим функцию от одного аргумента - списка.
Опять используем REPL для экпериментов с полученной функцией получения k-меров:
*BruteforceMotifSearch> k_mers 2 "ACTG" ["AC","CT","TG"] *BruteforceMotifSearch> k_mers 3 "AACCTTGG" ["AAC","ACC","CCT","CTT","TTG","TGG"]
Читатели, не имевшие опыта программирования на таких языках как Haskell наверняка скажут взглянув на код функции k_mers: “Это же ужасно не эффективно: сначала мы получем список ВСЕХ хвостов, часть из которых вовсе не нужна, потом формируем список из n-k+1 элементов, затем перебираем элементы уже этого списка, не лучше ли это сделать в одном цикле, вместо того, чтобы каждый раз передавать новый список на вход функциям?”. Это было бы справедливо, если бы не одно НО. Язык Haskell использует ленивый порядок вычислений. Это означает, что значения не вычисляются, пока они не понадобятся. Например, если мы взяли n-k+1 элемент списка и работаем с ними, то остальные элементы и не будут вычислены. Соответственно и список не формируется при вызове каждой из функций, а передается “обещание” вычислить этот список. А само вычисление происходит только в месте, где требуется получить конкретное значение. Про ленивые вычисления можно прочитать, например, здесь.
Теперь можно считать общее расстояние Хэминга (total Hamming distance):
total_dH seqs k_mer = sum $ map (min_dH k_mer) seqs where min_dH k_mer seq = minimum $ map (dH k_mer) $ k_mers (length k_mer) seq
Обратите внимание, что код функции total_dH отлично читается: суммировать минимальные расстояния для заданного k-мера для всех последовательностей.
Теперь осталось только сгенерировать список всех вариантов k-меров для перебора. Для этого воспользуемся функцией replicateM из модуля Control.Monad.
Допустим, нам необходимо получить все возможные 3-меры над алфавитом ДНК:
*BruteforceMotifSearch> replicateM 3 "ACGT" ["AAA","AAC","AAG","AAT","ACA","ACC","ACG","ACT","AGA","AGC","AGG","AGT","ATA","ATC","ATG","ATT","CAA","CAC","CAG","CAT","CCA","CCC","CCG","CCT","CGA","CGC","CGG","CGT","CTA","CTC","CTG","CTT","GAA","GAC","GAG","GAT","GCA","GCC","GCG","GCT","GGA","GGC","GGG","GGT","GTA","GTC","GTG","GTT","TAA","TAC","TAG","TAT","TCA","TCC","TCG","TCT","TGA","TGC","TGG","TGT","TTA","TTC","TTG","TTT"]
Собрав все описанные функции воедино и обозначив последовательность “ACGT” через dnaAlphabet, получим:
bruteforceMedianString seqs k = map snd $ filter minDist $ zip totalDistances allWords where totalDistances = map (total_dH seqs) allWords allWords = replicateM k dnaAlphabet minDist (n, _) = n == minimum totalDistances
Разберем эту функцию с “конца”:
для того, чтобы связать k-меры и их расстояния воспользуемся функцией zip, которая формирует из двух списков один список пар:
*BruteforceMotifSearch> zip ["AA", "AC", "AT"] [1,2,3] [("AA",1),("AC",2),("AT",3)]
Далее фильтруем полученный список находя элементы, которые равны минимальному общему расстоянию Хэмминга. Можно было бы обойтись просто поиском минимального расстояния и соответствующего k-мера, но в так мы не получим всех результатов с одинаковым минимальным значением. С увеличением k это будет все менее актуальным, т.к. значений с одинаковым расстоянием будет все меньше и меньше.
Для демонстрации работы воспользуемся случайно сгенерированными последовательностями, в которые “имплантирован” в случайной позиции каждой последовательности заданный k-мер.
Тесты и генерирование последовательностей в данной статье рассматриваться не будут, но с исходниками можно ознакомится здесь.
Генерируем 4 тестовые последовательности c “имплантированным” 7-мером AAACCCT:
CATCACAAAACCCTCGGTTACG AAAACCCTAGAACGGCCTGCGG ACAAGAAACCCTCAAGGACACT ACTAAACCCTCTGCTGAAGGTG
*MotifTests> bruteforceMedianString ["CATCACAAAACCCTCGGTTACG","AAAACCCTAGAACGGCCTGCGG","ACAAGAAACCCTCAAGGACACT","ACTAAACCCTCTGCTGAAGGTG"] 7 ["AAACCCT"]
Можно заметить, что поиск даже такого короткого k-мера затрачивается порядка 25 сек (если запускать в REPL на моем ноутбуке с Intel Core 2 Duo). Безусловно, если откомпилировать и включить опции оптимизации, то работать будет значительно быстрее, но тем не менее вычислительная сложность алгоритма характеризуется как O(4kkn), т.е. растет экспоненциально. Алгоритмы с такой сложностью по времени как правило не применимы на практике, но их можно улучшить с помощью уменьшения количества вариантов решений. Об этом я надеюсь рассказать в следующих статьях.
Список источников:
| < Предыдущая |
|---|
