Обновлено 05.03.2013
В продолжение предыдущей статьи рассмотрим решение задачи поиска регуляторных мотивов с помощью branch-and-bound алогоритма (на русском языке известен как “Метод ветвей и границ”), который позволяет существенно сократить время поиска решения.
Суть этого метода состоит в том, что в общем множестве решений (в котором осуществляется поиск методом полного перебора) отсекаются подмножества заведомо не оптимальных решений.
Для того, чтобы применить данный метод к задаче, необходимо проделать следующие шаги:
В прошлый раз, при использовании метода полного перебора, использовалась простая функция, которая генерировала все возможные варианты. В данном случае придется немного усложнить задачу и построить дерево решений.
Дерево должно выглядеть следующим образом:
Последний уровень дерева будет представлять собой все возможные варианты k-мера, высота дерева в этом случае будет k+1.
Для представления решений в виде дерева воспользуемся модулем Data.Tree. Чтобы сгенерировать дерево, воспользуемся функцией unfoldTree. Она принимает в качестве параметров функцию, которая генерирует узлы следующего уровня для заданного узла и значение корневого узла.
Функция, генерирующая узлы, будет выглядеть следующим образом:
nextLevel k prefix | length prefix < k = (prefix, map (appendChar prefix) dnaAlphabet) | otherwise = (prefix, []) where appendChar str c = str ++ [c]
Если длина переданного префикса меньше k, генерируем следующий уровень, в противоположном случае возвращаем признак окончания дерева.
Проверим, как это будет работать вместе с функцией unfoldTree. Для вывода дерева в пригодном для просмотра виде воспользуемся функцей drawTree, результат выполнения которой передадим putStr, которая и выведет на экран строку, представляющую дерево (результат немного изменен, чтобы занимать меньше места).
*BBMotifSearch> putStr $ drawTree $ unfoldTree (nextLevel 2) "" |
+- A
| +- AA
| +- AC
| +- AG
| `- AT
|
+- C
| +- CA
| +- CC
| +- CG
| `- CT
|
+- G
| +- GA
| +- GC
| +- GG
| `- GT
|
`- T
+- TA
+- TC
+- TG
`- TT
Рис. 2. Дерево для k=2
Кроме значений, представляющих собой строки длиной от 1 до k, в узлах дерева необходимо хранить также общее расстояние Хэмминга для каждой строки. Это позволит в дальнейшем применять процедуру отсечения. Для этого создадим специальный тип данных:
data TreeItem = Item { k_mer :: String , distance :: Int } deriving (Show, Eq)
Также определим процедуру сравнения для значений данного типа
instance Ord TreeItem where compare i1 i2 = compare (distance i1) (distance i2)
Будем сравнивать значения по их расстояниям.
Теперь можно определить функцию, возвращающую дерево решений:
searchTree seqs k = fmap f $ unfoldTree (nextLevel k) "" where f x = Item x (total_dH seqs x)
С работой unfoldTree мы уже знакомы, а вот функция fmap еще не встречалась. Эта функция декларирована для класса типов Functor. Она позволяет применить некоторую функцию к значению, содержащемуся в некотором контейнере, сам контейнер при этом не изменяется. Для списков это эквивалентно функции map. В случае деревьев переданная функция будет применяться к каждому узлу.
Таким образом, функция searchTree строит дерево из значений типа TreeItem, принимая в качестве параметров список исходных последовательностей, в которых ищется k-мер и длина k-мера.
Проверим в REPL результат работы функции. Для этого добавим определение последовательностей, которые мы будем использовать для примера:
sampleSeqs = ["TGACCGTGCCCTTGGA", "CCCTTGGAAGAAAAATGG", "AAACCTTGGACATGACT"]
Выражение, используемое для вывода дерева на экран, немного усложнится:
putStr $ drawTree $ fmap show $ searchTree sampleSeqs 2
Перед передачей результатов в функцию drawTree необходимо преобразовать узлы дерева в строковое представление, для чего и используется fmap. Функция show, примененная к значению типа TreeItem, даст его строковое представление.
Теперь получим результат в следующем виде:
Item {k_mer = "", distance = 0}
|
+- Item {k_mer = "A", distance = 0}
| +- Item {k_mer = "AA", distance = 1}
| +- Item {k_mer = "AC", distance = 1}
| +- Item {k_mer = "AG", distance = 2}
| `- Item {k_mer = "AT", distance = 1}
|
+- Item {k_mer = "C", distance = 0}
| +- Item {k_mer = "CA", distance = 2}
| +- Item {k_mer = "CC", distance = 0}
| +- Item {k_mer = "CG", distance = 2}
| `- Item {k_mer = "CT", distance = 0}
|
+- Item {k_mer = "G", distance = 0}
| +- Item {k_mer = "GA", distance = 0}
| +- Item {k_mer = "GC", distance = 2}
| +- Item {k_mer = "GG", distance = 0}
| `- Item {k_mer = "GT", distance = 2}
|
`- Item {k_mer = "T", distance = 0}
+- Item {k_mer = "TA", distance = 3}
+- Item {k_mer = "TC", distance = 3}
+- Item {k_mer = "TG", distance = 0}
`- Item {k_mer = "TT", distance = 0}
Рис. 3. Дерево хранящее в узлах значения типа TreeItem(k=2)
Теперь можно приступать к формированию оценки “перспективности” ветвей дерева. Подробнее об этом можно прочитать в [2], здесь мы воспользуемся результатами, изложенными в книге.
Критерий, по которому определяется необходимость двигаться по ветке дерева, прост:
если расстояние в узле дерева больше некоторого текущего расстояния, то двигаться ниже не имеет смысла, т.к. узлы, лежащие ниже, улучшить расстояние не смогут. В них может быть либо то же самое расстояние, либо большее, что нам не подходит, т.к. решается задача минимизации.
Значит, осталось только двигаться по дереву и разворачивать ветки, руководствуясь этим критерием. Обход дерева должен осуществляться в прямом порядке (preorder), он хорошо проиллюстрирован в [2], а также в презентации, так что на нем подробно останавливаться не будем.
Чтобы реализовать этот алгоритм необходимы некоторые вспомогательные функции для работы с деревом.
В работе [3] приведены некоторые функции для работы с деревом, которые будут полезны для реализации branch-and-bound алгоритма. После небольшой их модификации и добавления некоторых функций получился модуль TreeUtils, в котором и собраны все вспомогательные функции, которые нам понадобятся.
Основные из них:
levels' :: Tree a -> [Forest a]
работает почти как оригинальная функция levels из модуля Data.Tree, за исключением того, что возвращает не списки значений(roolLabel) узлов дерева по уровням, а список значений типа Forest (Forest a эквивалетно [Tree a]).
Например(оператор !! используется для получения элемента по индексу),
(levels’ $ searchTree sampleSeqs 2)!!1
Вернет список с содержимым 1-го уровня (если считать корень нулевым уровнем):
[Item {k_mer = "A", distance = 0},Item {k_mer = "C", distance = 0},Item {k_mer = "G", distance = 0},Item {k_mer = "T", distance = 0}]
Следующая важная функция - это функция prune. С её помощью можно фильтровать дерево, отрезая ветки, которые не содержат оптимальных решений.
Поэкспериментируем с функцией prune в REPL. Например, обрежем ветки в узлах с расстоянием ≥ 2.
*BBMotifSearch> let prunedTree = prune (\x -> distance x >=2) $ searchTree sampleSeqs 2
*BBMotifSearch> putStr $ drawTree$ fmap show $ prunedTree
Item {k_mer = "", distance = 0}
|
+- Item {k_mer = "A", distance = 0}
| |
| +- Item {k_mer = "AA", distance = 1}
| |
| +- Item {k_mer = "AC", distance = 1}
| |
| `- Item {k_mer = "AT", distance = 1}
|
+- Item {k_mer = "C", distance = 0}
| |
| +- Item {k_mer = "CC", distance = 0}
| |
| `- Item {k_mer = "CT", distance = 0}
|
+- Item {k_mer = "G", distance = 0}
| |
| +- Item {k_mer = "GA", distance = 0}
| |
| `- Item {k_mer = "GG", distance = 0}
|
`- Item {k_mer = "T", distance = 0}
|
+- Item {k_mer = "TG", distance = 0}
|
`- Item {k_mer = "TT", distance = 0}
Рис. 4. Дерево, у которого отсечены ветки с расстоянием ≥2
Если сравнить полученный результат с исходным деревом (Рис.3.), можно увидеть сократившееся количество ветвей.
Теперь необходимо организовать обход дерева с организацией отсечения ненужных ветвей. Мы немного “схитрим”, воспользовавшись для дальнейших манипуляций с деревом функцией levels’ вместо явного обхода дерева в прямом порядке (preorder).
Введем еще одну функцию, которая будет возвращать вершины нижнего уровня (листья), которые еще доступны для обхода.
availableVertices tree k
| null lastLevels || null nextToLastVertices = []
| otherwise = subForest $ head nextToLastVertices
where
lastLevels = drop (k-1) (levels' tree)
nextToLastVertices = dropWhile (null . subForest) $ head lastLevels
Механизм работы функции таков:
Теперь можно легко определить, что поиск окончен: availableVertices вернет пустой список, т.е не осталось вершин для дальнейшего рассмотрения.
Теперь скомпонуем финальную функцию, которая будет находить потенциальное решение, фильтровать дерево решений и рекурсивно продолжать поиск до тех пор, пока не останется вершин, которые могут содержать лучшее решение.
improve k sTree | null $ availableVertices prunedTree k = k_mer minItem | otherwise = improve k prunedTree where minItem = minimum $ map rootLabel nextVertices prunedTree = prune (>= minItem) sTree nextVertices = availableVertices sTree k
Естественно, что такой способ решения задачи работает только для "ленивых" языков, иначе бы все дерево решений вычислялось. В случае Haskell имеем как раз ленивую модель вычислений, следовательно дерево будет строится по необходимости и, если условие разворачивания ветвей не выполнилось, то вершины лежащие ниже просто не будут вычисляться.
Для тестирования воспользуемся теми же данными, что и для bruteforce - алгоритма из предыдущей статьи(ссылка):
*BBMotifSearch> let sampleSeqs = ["CATCACAAAACCCTCGGTTACG","AAAACCCTAGAACGGCCTGCGG","ACAAGAAACCCTCAAGGACACT","ACTAAACCCTCTGCTGAAGGTG"]
*BBMotifSearch> bbMedianString sampleSeqs 7
"AAACCCT"
Результат совпадает с данными из предыдущей статьи. Единственным отличием является то, что решение представлено одним k-мером(первым из оптимальных), тогда как bruteforce - алгоритм находит все оптимальные решения. Но это различие, как было отмечено в предыдущей статье, перестает быть актуальным с увеличением k.
Теперь посмотрим, что же улучшилось в алгоритме и как это повлияло на его скорость. Для этого воспользуемся профайлером Haskell. Чтобы скомпилировать исходный код с возможностью профилирования необходимо указать соответствующие опции компилятора, например:
ghc -rtsopts -prof -fprof-auto -O2 -main-is BBMotifSearch BBMotifSearch.hs
Далее запускаем программу, указав ей необходимые параметры и опции для запуска профайлера:
./BBMotifSearch purT_orthologs.fas 8 bb +RTS -p
Здесь purT_orthologs.fas - файл с последовательностями в формате FASTA, 8 - длина мотива для поиска, bb - алгоритм (bb - branch-and-bound, bf - bruteforce) +RTS - блок опций специальных параметров, -p - включение профайлера (после завершения программы результат профилирования будет записан в файл BBMotifSearch.prof). В качестве исходных данных будем использовать 55 последовательностей длиной 200 символов каждая. Исходные данные также можно найти в репозитории (файл bio/MotifSearch/purT_orthologs.fas).
Для чтения формата FASTA используется модули Bio.Sequence.Fasta, Bio.Sequence.SeqData из библиотеки bio. Также в код была добавлена функция main, которая должна обязательно присутствовать в коде , если планируется осуществлять компиляцию в исполнимый файл. Именно она будет просматривать параметры командной строки и читать файл с данными и осуществлять вывод результата. Т.е. все побочные эффекты отделены от остального кода, который представляет собой чистые функции(подробнее, о том как должна быть структурирована программа на Haskell можно ознакомится, например, в [1]).
Ищем 8-мер bruteforce - алгоритмом
./BBMotifSearch purT_orthologs.fas 8 bf +RTS -p
total time = 912.31 secs (912310 ticks @ 1000 us, 1 processor)
total alloc = 570,633,453,300 bytes (excludes profiling overheads)
COST CENTRE MODULE %time %alloc
dH BruteforceMotifSearch 54.5 58.5
k_mers BruteforceMotifSearch 23.5 36.1
dH.diff BruteforceMotifSearch 13.7 0.0
total_dH.min_dH BruteforceMotifSearch 6.6 5.4
k_mers.n BruteforceMotifSearch 1.5 0.0
individual inherited
COST CENTRE MODULE no. entries %time %alloc %time %alloc
MAIN MAIN 59 0 0.0 0.0 100.0 100.0
main BBMotifSearch 119 0 0.0 0.0 100.0 100.0
seqdata Bio.Sequence.SeqData 137 55 0.0 0.0 0.0 0.0
bruteforceMedianString BruteforceMotifSearch 121 1 0.0 0.0 100.0 100.0
bruteforceMedianString.minDist BruteforceMotifSearch 127 65536 0.0 0.0 0.0 0.0
bruteforceMedianString.totalDistances BruteforceMotifSearch 126 1 0.0 0.0 100.0 100.0
total_dH BruteforceMotifSearch 128 65536 0.1 0.0 100.0 100.0
total_dH.min_dH BruteforceMotifSearch 133 3604480 6.6 5.4 99.9 100.0
dH BruteforceMotifSearch 142 695664640 54.5 58.5 68.3 58.5
dH.diff BruteforceMotifSearch 143 5565317120 13.7 0.0 13.7 0.0
k_mers BruteforceMotifSearch 134 3604480 23.5 36.1 25.0 36.1
k_mers.n BruteforceMotifSearch 135 3604480 1.5 0.0 1.5 0.0
bruteforceMedianString.words BruteforceMotifSearch 122 1 0.0 0.0 0.0 0.0
readFasta Bio.Sequence.Fasta 120 1 0.0 0.0 0.0 0.0
mkSeqs Bio.Sequence.Fasta 130 0 0.0 0.0 0.0 0.0
mkSeq Bio.Sequence.Fasta 138 55 0.0 0.0 0.0 0.0
mkSeq.isSeq Bio.Sequence.Fasta 139 55 0.0 0.0 0.0 0.0
blocks Bio.Sequence.Fasta 132 0 0.0 0.0 0.0 0.0
CAF GHC.IO.Encoding.Iconv 99 0 0.0 0.0 0.0 0.0
CAF Text.Read.Lex 97 0 0.0 0.0 0.0 0.0
CAF GHC.IO.FD 86 0 0.0 0.0 0.0 0.0
CAF GHC.IO.Encoding 84 0 0.0 0.0 0.0 0.0
CAF GHC.Conc.Signal 80 0 0.0 0.0 0.0 0.0
CAF GHC.IO.Handle.FD 76 0 0.0 0.0 0.0 0.0
CAF Bio.Sequence.Fasta 70 0 0.0 0.0 0.0 0.0
mkSeq Bio.Sequence.Fasta 140 0 0.0 0.0 0.0 0.0
mkSeq.isSeq Bio.Sequence.Fasta 141 0 0.0 0.0 0.0 0.0
blocks Bio.Sequence.Fasta 131 1 0.0 0.0 0.0 0.0
mkSeqs Bio.Sequence.Fasta 129 1 0.0 0.0 0.0 0.0
CAF Bio.Sequence.SeqData 69 0 0.0 0.0 0.0 0.0
toStr Bio.Sequence.SeqData 136 1 0.0 0.0 0.0 0.0
CAF BruteforceMotifSearch 68 0 0.0 0.0 0.0 0.0
bruteforceMedianString BruteforceMotifSearch 124 0 0.0 0.0 0.0 0.0
bruteforceMedianString.words BruteforceMotifSearch 125 0 0.0 0.0 0.0 0.0
dnaAlphabet BruteforceMotifSearch 123 1 0.0 0.0 0.0 0.0
CAF BBMotifSearch 66 0 0.0 0.0 0.0 0.0
main BBMotifSearch 118 1 0.0 0.0 0.0 0.0
main BBMotifSearch 144 0 0.0 0.0 0.0 0.0
readFasta Bio.Sequence.Fasta 145 0 0.0 0.0 0.0 0.0
Ищем 8-мер branch-and-bound алгоритмом
./BBMotifSearch purT_orthologs.fas 8 bb +RTS -p
total time = 340.82 secs (340820 ticks @ 1000 us, 1 processor)
total alloc = 170,229,832,956 bytes (excludes profiling overheads)
COST CENTRE MODULE %time %alloc
dH BruteforceMotifSearch 56.0 61.1
dH.diff BruteforceMotifSearch 21.6 1.6
k_mers BruteforceMotifSearch 14.3 27.8
total_dH.min_dH BruteforceMotifSearch 7.1 9.5
individual inherited
COST CENTRE MODULE no. entries %time %alloc %time %alloc
MAIN MAIN 59 0 0.0 0.0 100.0 100.0
main BBMotifSearch 119 0 0.0 0.0 100.0 100.0
seqdata Bio.Sequence.SeqData 147 55 0.0 0.0 0.0 0.0
bbMedianString BBMotifSearch 121 1 0.0 0.0 100.0 100.0
bbMedianString.sTree BBMotifSearch 130 1 0.0 0.0 100.0 100.0
searchTree BBMotifSearch 131 1 0.0 0.0 100.0 100.0
searchTree.f BBMotifSearch 136 17168 0.0 0.0 100.0 100.0
total_dH BruteforceMotifSearch 137 17168 0.1 0.0 100.0 100.0
total_dH.min_dH BruteforceMotifSearch 142 944240 7.1 9.5 99.9 99.9
dH BruteforceMotifSearch 152 183525760 56.0 61.1 77.5 62.7
dH.diff BruteforceMotifSearch 153 1217379240 21.6 1.6 21.6 1.6
k_mers BruteforceMotifSearch 144 944240 14.3 27.8 15.2 27.8
k_mers.n BruteforceMotifSearch 145 944240 0.9 0.0 0.9 0.0
nextLevel BBMotifSearch 133 17169 0.0 0.0 0.0 0.0
nextLevel.appendChar BBMotifSearch 143 17168 0.0 0.0 0.0 0.0
improve BBMotifSearch 122 8 0.0 0.0 0.0 0.0
improve.minItem BBMotifSearch 159 8 0.0 0.0 0.0 0.0
compare BBMotifSearch 160 3 0.0 0.0 0.0 0.0
improve.nextVertices BBMotifSearch 154 8 0.0 0.0 0.0 0.0
availableVertices BBMotifSearch 155 8 0.0 0.0 0.0 0.0
availableVertices.nextToLastVertices BBMotifSearch 158 8 0.0 0.0 0.0 0.0
availableVertices.lastLevels BBMotifSearch 156 8 0.0 0.0 0.0 0.0
levels' TreeUtils 157 8 0.0 0.0 0.0 0.0
improve.prunedTree BBMotifSearch 126 8 0.0 0.0 0.0 0.0
prune TreeUtils 127 8 0.0 0.0 0.0 0.0
filterTree TreeUtils 128 8 0.0 0.0 0.0 0.0
foldTree TreeUtils 129 112963 0.0 0.0 0.0 0.0
filterTree.f TreeUtils 132 112963 0.0 0.0 0.0 0.0
compare BBMotifSearch 135 112955 0.0 0.0 0.0 0.0
availableVertices BBMotifSearch 123 8 0.0 0.0 0.0 0.0
availableVertices.nextToLastVertices BBMotifSearch 163 8 0.0 0.0 0.0 0.0
availableVertices.lastLevels BBMotifSearch 124 8 0.0 0.0 0.0 0.0
levels' TreeUtils 125 8 0.0 0.0 0.0 0.0
readFasta Bio.Sequence.Fasta 120 1 0.0 0.0 0.0 0.0
mkSeqs Bio.Sequence.Fasta 139 0 0.0 0.0 0.0 0.0
mkSeq Bio.Sequence.Fasta 148 55 0.0 0.0 0.0 0.0
mkSeq.isSeq Bio.Sequence.Fasta 149 55 0.0 0.0 0.0 0.0
blocks Bio.Sequence.Fasta 141 0 0.0 0.0 0.0 0.0
CAF GHC.IO.Encoding.Iconv 99 0 0.0 0.0 0.0 0.0
CAF Text.Read.Lex 97 0 0.0 0.0 0.0 0.0
CAF GHC.IO.FD 86 0 0.0 0.0 0.0 0.0
CAF GHC.IO.Encoding 84 0 0.0 0.0 0.0 0.0
CAF GHC.Conc.Signal 80 0 0.0 0.0 0.0 0.0
CAF GHC.IO.Handle.FD 76 0 0.0 0.0 0.0 0.0
CAF Bio.Sequence.Fasta 70 0 0.0 0.0 0.0 0.0
mkSeq Bio.Sequence.Fasta 150 0 0.0 0.0 0.0 0.0
mkSeq.isSeq Bio.Sequence.Fasta 151 0 0.0 0.0 0.0 0.0
blocks Bio.Sequence.Fasta 140 1 0.0 0.0 0.0 0.0
mkSeqs Bio.Sequence.Fasta 138 1 0.0 0.0 0.0 0.0
CAF Bio.Sequence.SeqData 69 0 0.0 0.0 0.0 0.0
toStr Bio.Sequence.SeqData 146 1 0.0 0.0 0.0 0.0
CAF BruteforceMotifSearch 68 0 0.0 0.0 0.0 0.0
dnaAlphabet BruteforceMotifSearch 134 1 0.0 0.0 0.0 0.0
CAF BBMotifSearch 66 0 0.0 0.0 0.0 0.0
main BBMotifSearch 118 1 0.0 0.0 0.0 0.0
main BBMotifSearch 161 0 0.0 0.0 0.0 0.0
readFasta Bio.Sequence.Fasta 162 0 0.0 0.0 0.0 0.0
При k=8 разница во времени работы составит: tbruteforce/tb&b=912.31/340.82≈2.7. Эта разница будет расти с увеличением k.
В выдаче профайлера также можно почерпнуть и другую полезную информацию - количество вызовов различных функций, их вклад в общее время выполнения, а также данные по распределению памяти.
Сравним количество вызовов основных функций:
| Функция | bruteforce | branch-and-bound |
|---|---|---|
| total_dH |
65 536 |
17 168 |
| dH |
695 664 640 |
183 525 760 |
| k_mers | 3 604 480 | 944 240 |
Видно, что количество вызовов во втором случае значительно меньше. По количеству вызовов total_dH можно судить о том, сколько узлов было просмотрено у дерева решений.
Данное решение обладает некоторой универсальностью, позволяющей изменить branch-and-bound стратегию без изменения самого поискового механизма. Для этого достаточно будет поменять механизм формирования дерева решений.
Здесь, например, приведен пример использования более “агрессивной” стратегии, позволяющей отсечь больше ветвей, не содержащих оптимальных решений.
Реализацию этой стратегии можно посмотреть в репозитории исходного кода (файл bio/MotifSearch/BBMotifSearch1.hs).
Также возможно получить значительное ускорение работы(порядка 6 раз для k=8) используя тип ByteString вместо стандартного типа String. Результат также представлен в репозитории (файл bio/MotifSearch/BBMotifSearchByteString.hs).
UPD1: исправлена выдача профайлера для bruteforce-алгоритма. При использовании GHC 7.0.3 неправильно считалось количество вызовов функций. После компиляции GHC 7.4.1 был получен правильный результат. К аргументам командной строки добавлен параметр, задающий алгоритм: bb - branch-and-bound, bf - bruteforce.
UPD2: в исходный код добавлена реализация обхода дерева в прямом (preorder) порядке с "протаскиванием" состояния (минимального, на данный момент, общего расстояния Хэмминга) с помощью State Monad. Монады в Haskell требуют отдельных разъяснений, которые выходят за рамки этой статьи. Про монады можно почитать, например, здесь, здесь и здесь. Также добавлено значение bb-preorder для задания этого способа обхода через командную строку.
| Следующая > |
|---|