Решение задачи о поиске наиболее часто встречающейся подпоследовательности с заданным максимальным числом мутаций на haskell

Эта задача имеет большое значение в биологии и часто формулируется как поиск наиболее часто встречающегося k-мера с максимально допустимым числом мутаций n внутри заданной генетической последовательности. Генетическая последовательность обычно моделируется строкой произвольной длины, составленной из символов некоторого алфавита. Соответственно, k-мер — это подстрока длиной k внутри генетической последовательности. В постановке задачи искомые k-меры могут накладываться друг на друга. Так, если общая длина исходной последовательности равна 10, то в ней определены один 10-мер, два 9-мера, три 8-мера и т.д. Благодаря условию о допустимых мутациях искомая наиболее часто повторяющаяся подпоследовательность может вообще ни разу не встретиться в исходной последовательности! Для решения задачи предложим такой алгоритм: для каждого k-мера в заданной последовательности находим все его мутации, включая исходное значение, складываем их вместе и находим наиболее часто встречающийся элемент: он и будет ответом на поставленную задачу. Я приведу решение, а ниже его прокомментирую.

import qualified Data.List as L
import Data.Ord (comparing)
import Control.Arrow

alphabet :: [Char]
alphabet = "ACGT"

subSeqs :: Int -> String -> [String]
subSeqs n = takeWhile ((==n) . length) . map (take n) . L.tails

nearSeqs :: Int -> String -> [String]
nearSeqs n = map head . L.group . L.sort . nearSeqs' n
    where nearSeqs' 0 s = [s]
          nearSeqs' n s =
            concatMap (\(a, b) -> map (a++) b) $
                concatMap (\m -> [(z ++ [x], nearSeqs' (n - 1) $ tail z') |
                                  x <- alphabet, let (z, z') = splitAt m s])
                [0 .. length s - 1]

grpSeqs :: Int -> [String] -> [[String]]
grpSeqs n = L.sortBy (flip $ comparing length) . L.group . L.sort . seqs
    where seqs = concatMap $ nearSeqs n

main = do
    let seqs = ["ACGTTGCATGTCGCATGATGCATGAGAGCT",
                "CACAGTAGGCGCCGGCACACACAGCCCCGGGCCCCGGGCCGCCCCGGGCCGGC\
                \GGCCGCCGGCGCCGGCACACCGGCACAGCCGTACCGGCACAGTAGTACCGGCC\
                \GGCCGGCACACCGGCACACCGGGTACACACCGGGGCGCACACACAGGCGGGCG\
                \CCGGGCCCCGGGCCGTACCGGGCCGCCGGCGGCCCACAGGCGCCGGCACAGTA\
                \CCGGCACACACAGTAGCCCACACACAGGCGGGCGGTAGCCGGCGCACACACAC\
                \ACAGTAGGCGCACAGCCGCCCACACACACCGGCCGGCCGGCACAGGCGGGCGG\
                \GCGCACACACACCGGCACAGTAGTAGGCGGCCGGCGCACAGCC"]
    {-print $ nearSeqs 2 "TTTT"-}
    {-mapM_ print' $ take 10 $ grpSeqs 1 $ subSeqs 4 $ head seqs-}
    mapM_ print' $ take 10 $ grpSeqs 2 $ subSeqs 10 $ seqs !! 1
    where print' = print . (head &&& length)

В функции main происходит тестирование алгоритма. Функция alphabet задает алфавит последовательности. Как видим, он состоит всего из четырех символов. Я намеренно задал тип alphabet как [Char], а не String, дабы подчеркнуть, что это алфавит, а не простая строка. Для простоты я не стал каким-либо образом верифицировать исходные последовательности в seqs, объявленной в main. Функция subSeqs возвращает все k-меры из заданной последовательности (ее первый формальный параметр n соответствует числу k). Функция nearSeqs возвращает список уникальных мутаций (включая исходное значение) некоторой последовательности, ее формальный параметр n задает максимально допустимое количество мутаций. Реализация этой функции основана на рекурсивном применении генератора списков по всевозможным перестановкам общего числа n позиций внутри заданной последовательности. Для каждой перестановки генерируются всевозможные мутации на основе алфавита alphabet: при таком подходе результирующий список может содержать повторяющиеся элементы, поэтому он фильтруется функцией map head . L.group . L.sort. Как выяснилось впоследствии, правильный выбор этой функции очень сильно влияет на скорость программы. Так, в моем оригинальном варианте L.nub . L.sort программа работала в два раза медленнее! Просто L.nub работает еще медленнее (все-таки это функция класса O(n^2)). Текущий вариант я подсмотрел здесь. Позже я нашел самый быстрый вариант: Data.List.Ordered.nubSort, но все же решил на нем не останавливаться, поскольку он требует импорта еще одного модуля. К слову, функцию nearSeqs можно определить без необходимости финальной сортировки списка, заранее вынимая те элементы из alpahabet, которые приведут к дубликатам.

nearSeqs'' :: Int -> String -> [String]
nearSeqs'' n s = s : concatMap (`nearSeqs'` s) [1 .. n]
    where nearSeqs' 0 s = [s]
          nearSeqs' n s =
            concatMap (\(a, b) -> map (a++) b) $
                concatMap (\m -> [(z ++ [x], nearSeqs' (n - 1) $ tail z') |
                                  let (z, z') = splitAt m s,
                                  x <- L.delete (head z') alphabet])
                [0 .. length s - 1]

Однако, эта функция работает медленнее, поэтому я ее забраковал. Чтобы посмотреть на функцию nearSeqs в деле, раскомментируйте первую закомментированную строку в main, скомпилируйте программу и запустите ее. Функция grpSeqs вызывает nearSeqs, сортирует и группирует возвращенный ею список, и далее сортирует сгруппированный список по длине его подсписков в порядке убывания. В функции main из этого списка берется 10 первых подсписков и выводятся их первый элемент (все элементы подсписков должны совпадать) и длина. Для компиляции всех примеров в этой статье я использовал флаг -O2, имя исходника — gen.hs. Вывод программы с замером времени выполнения:

time ./gen
("GCACACAGAC",20)
("GCGCACACAC",20)
("ACACACACAC",19)
("CCCGCACACA",19)
("CACACACACG",18)
("CGCGCACACA",18)
("CGGCACACAC",18)
("GCACACACCC",18)
("GCACACACGC",18)
("GCCCGCACAC",18)

real    0m0.639s
user    0m0.612s
sys     0m0.022s

Чуть больше, чем полсекунды. Неплохо. Однако, этого я достиг не сразу и поэтому перепробовал множество подходов. В частности, я пробовал использовать для сборки результатов nearSeqs хэш-таблицу Data.HashTable.IO. В список импортируемых модулей добавляем

import qualified Data.HashTable.IO as H
import Data.Maybe (fromMaybe)

В определения записываем

type HashTable k v = H.CuckooHashTable k v

grpSeqsHt :: Int -> Int -> [String] -> IO (HashTable String Int)
grpSeqsHt m n s = do
    h <- H.newSized m
    let insert' h k = H.lookup h k >>= H.insert h k . (+1) . fromMaybe 0
    mapM_ (mapM_ (insert' h) . nearSeqs n) s
    return h

mostOftenSeqsHt :: HashTable String Int -> IO [(String, Int)]
mostOftenSeqsHt = H.foldM findMax []
    where findMax [] a = return [a]
          findMax c@((_, v') : _) a@(_, v) =
            return $ case v `compare` v' of GT -> [a]; EQ -> a : c; LT -> c

В функцию main записываем строку

    grpSeqsHt 70000 2 (subSeqs 10 $ seqs !! 1) >>= mostOftenSeqsHt >>= print

Функция grpSeqsHt создает хэш-таблицу исходного размера m (величина 70000 подобрана эмпирическим путем: примерно столько уникальных k-меров генерируется из seqs !! 1). Ключами таблицы являются значения k-меров, значениями — их общее количество в исходных данных. Если при добавлении новой ячейки ее еще не существует в таблице (это проверяется функцией H.lookup), то она создается со значением 1, иначе соответствующее ей значение увеличивается на единицу. Функция mostOftenSeqsHt проходится по уже заполненной таблице и возвращает список наиболее часто встречающихся k-меров (заметьте отличие: первый вариант решения задачи формировал список всех k-меров в порядке уменьшения их частоты в исходной последовательности). Второй вариант с использованием хэш-таблицы учитывает тот факт, что вычисление списка наиболее часто встречающихся k-меров можно произвести прямо во время заполнения таблицы внутри функции grpSeqsHt, и функция mostOftenSeqsHt становится ненужной. Итак, добавляем в список импортируемых модулей

import Control.Monad (foldM)

, в определения

grpSeqsHt' :: Int -> Int -> [String] ->
    IO (HashTable String Int, [(String, Int)])
grpSeqsHt' m n s = do
    h <- H.newSized m
    let insert' h k =
            H.lookup h k >>=
                (\v -> H.insert h k v >> return (k, v)) . (+1) . fromMaybe 0
        findMax [] a = return [a]
        findMax c@((_, v') : _) a@(_, v) =
            return $ case v `compare` v' of GT -> [a]; EQ -> a : c; LT -> c
    c <- foldM (\c s -> foldM
                    (\c k -> insert' h k >>= findMax c) c $ nearSeqs n s) [] s
    return (h, c)

, и в функцию main

    grpSeqsHt' 70000 2 (subSeqs 10 $ seqs !! 1) >>= print . snd

В этом случае использование хэш-таблицы становится побочным артефактом реализации, поскольку за пределами функции grpSeqsHt’ она не используется, несмотря на то, что grpSeqsHt’ возвращает ее как первый элемент кортежа. Как ни странно, варианты с хэш-таблицей оказались медленнее. Так, первый вариант выдал

time ./gen 
[("GCACACAGAC",20),("GCGCACACAC",20)]

real    0m0.656s
user    0m0.643s
sys     0m0.011s

, а второй

time ./gen 
[("GCGCACACAC",20),("GCACACAGAC",20)]

real    0m0.676s
user    0m0.664s
sys     0m0.009s

Ну и, собственно, чемпион! Возьмем исходное решение и чуточку изменим алгоритм. Предположим, что список, возвращенный из функции subSeqs, может содержать дубликаты (как оказалось, seqs !! 1 действительно содержит 98 дубликатов k-меров при их общем числе 352). Для дубликатов не нужно получать список мутаций каждый раз заново: достаточно получить список единожды и реплицировать его по числу дубликатов. Давайте просто перепишем функцию grpSeqs в предположении, что она по-прежнему будет получать данные из subSeqs.

grpSeqs :: Int -> [String] -> [[String]]
grpSeqs n = L.sortBy (flip $ comparing length) . L.group . L.sort . seqs
    where seqs = concatMap (\(s, m) -> concat . replicate m . nearSeqs n $ s) .
                    map (head &&& length) . L.group . L.sort

Вывод программы:

time ./gen 
("GCACACAGAC",20)
("GCGCACACAC",20)
("ACACACACAC",19)
("CCCGCACACA",19)
("CACACACACG",18)
("CGCGCACACA",18)
("CGGCACACAC",18)
("GCACACACCC",18)
("GCACACACGC",18)
("GCCCGCACAC",18)

real    0m0.546s
user    0m0.525s
sys     0m0.018s