как сопоставить образец последовательности ДНК

Пояснение алгоритма:

• Имея последовательность S с символами s(1), s(2),…, s(N).
• Пусть B(i) — наилучшая сжатая последовательность с элементами s(1), s(2),…,s(i).
• Так, например, B(3) будет лучшей сжатой последовательностью для s(1), s(2), s(3).
• Нам нужно знать B(N).

Чтобы найти его, будем действовать по индукции. Мы хотим вычислить B(i+1), зная B(i), B(i-1), B(i-2), …, B(1), B(0), где B(0) пусто последовательность, и и B(1) = s(1). В то же время это является доказательством оптимальности решения. ;)

Чтобы вычислить B(i+1), мы выберем лучшую последовательность среди кандидатов:

1. Последовательности-кандидаты, в которых последний блок имеет один элемент:

   B(i)s(i+1)1 B(i-1)s(i+1)2 ; только если s(i) = s(i+1) B(i-2)s(i+1)3 ; только если s(i-1) = s(i) и s(i) = s(i+1) … B(1)s(i+1)[i-1] ; только если s(2)=s(3) и s(3)=s(4) и … и s(i) = s(i+1) B(0)s(i+1)i = s(i +1)я; только если s(1)=s(2) и s(2)=s(3) и … и s(i) = s(i+1)

2. Последовательности-кандидаты, в которых последний блок состоит из 2 элементов:

   B(i-1)s(i)s(i+1)1 B(i-3)s(i)s(i+1)2 ; только если s(i-2)s(i-1)=s(i)s(i+1) B(i-5)s(i)s(i+1)3 ; только если s(i-4)s(i-3)=s(i-2)s(i-1) и s(i-2)s(i-1)=s(i)s(i+1) ) …

3. Последовательности-кандидаты, в которых последний блок состоит из 3 элементов:

   …

4. Последовательности-кандидаты, в которых последний блок состоит из 4 элементов:

   …

   …

5. Последовательности-кандидаты, в которых последний блок содержит n+1 элементов:

   s(1)s(2)s(3)………s(i+1)

Для каждой возможности алгоритм останавливается, когда блок последовательности больше не повторяется. И это все.

Алгоритм будет примерно таким в коде psude-c:

B(0) = “”
for (i=1; i<=N; i++) {
    // Calculate all the candidates for B(i)
    BestCandidate=null
    for (j=1; j<=i; j++) {
        Calculate all the candidates of length (i)

        r=1;
        do {
             Candidadte = B([i-j]*r-1) s(i-j+1)…s(i-1)s(i) r
             If (   (BestCandidate==null) 
                      || (Candidate is shorter that BestCandidate)) 
                 {
            BestCandidate=Candidate.
                 }
             r++;
        } while (  ([i-j]*r <= i) 
             &&(s(i-j*r+1) s(i-j*r+2)…s(i-j*r+j) == s(i-j+1) s(i-j+2)…s(i-j+j))

    }
    B(i)=BestCandidate
}

Надеюсь, что это может помочь немного больше.

Полная программа на C, выполняющая требуемую задачу, приведена ниже. Он работает за O (n ^ 2). Центральная часть состоит всего из 30 строк кода.

РЕДАКТИРОВАНИЕ Я немного реструктурировал код, изменил имена переменных и добавил комментарий, чтобы сделать его более читабельным.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <limits.h>


// This struct represents a compressed segment like atg4, g3,  agc1
    struct Segment {
        char *elements;
        int nElements;
        int count;
    };
         // As an example, for the segment agagt3  elements would be:
         // {
         //      elements: "agagt",
         //      nElements: 5,
         //      count: 3
         // }

    struct Sequence {
        struct Segment lastSegment;
        struct Sequence *prev;      // Points to a sequence without the last segment or NULL if it is the first segment
        int totalLen;  // Total length of the compressed sequence.
    };
       // as an example, for the sequence agt32ta5, the representation will be:
       // {
       //     lastSegment:{"ta" , 2 , 5},
       //     prev: @A,
       //     totalLen: 8
       // }
       // and A will be
       // {
       //     lastSegment{ "agt", 3, 32},
       //     prev: NULL,
       //     totalLen: 5
       // }


// This function converts a sequence to a string.
// You have to free the string after using it.
// The strategy is to construct the string from right to left.

char *sequence2string(struct Sequence *S) {
    char *Res=malloc(S->totalLen + 1);
    char *digits="0123456789";

    int p= S->totalLen;
    Res[p]=0;

    while (S!=NULL) {
            // first we insert the count of the last element.
            // We do digit by digit starting with the units.
            int C = S->lastSegment.count;
            while (C) {
                p--;
                Res[p] = digits[ C % 10 ];
                C /= 10;
            }

            p -= S->lastSegment.nElements;
            strncpy(Res + p , S->lastSegment.elements, S->lastSegment.nElements);

            S = S ->prev;
    }


    return Res;
}


// Compresses a dna sequence.
// Returns a string with the in sequence compressed.
// The returned string must be freed after using it.
char *dnaCompress(char *in) {
    int i,j;

    int N = strlen(in);;            // Number of elements of a in sequence.



    // B is an array of N+1 sequences where B(i) is the best compressed sequence sequence of the first i characters.
    // What we want to return is B[N];
    struct Sequence *B;
    B = malloc((N+1) * sizeof (struct Sequence));

    // We first do an initialization for i=0

    B[0].lastSegment.elements="";
    B[0].lastSegment.nElements=0;
    B[0].lastSegment.count=0;
    B[0].prev = NULL;
    B[0].totalLen=0;

    // and set totalLen of all the sequences to a very HIGH VALUE in this case N*2 will be enougth,  We will try different sequences and keep the minimum one.
    for (i=1; i<=N; i++) B[i].totalLen = INT_MAX;   // A very high value

    for (i=1; i<=N; i++) {

        // at this point we want to calculate B[i] and we know B[i-1], B[i-2], .... ,B[0]
        for (j=1; j<=i; j++) {

            // Here we will check all the candidates where the last segment has j elements

            int r=1;                  // number of times the last segment is repeated
            int rNDigits=1;           // Number of digits of r
            int rNDigitsBound=10;     // We will increment r, so this value is when r will have an extra digit.
                                      // when r = 0,1,...,9 => rNDigitsBound = 10
                                      // when r = 10,11,...,99 => rNDigitsBound = 100
                                      // when r = 100,101,.,999 => rNDigitsBound = 1000 and so on.

            do {

                // Here we analitze a candidate B(i).
                // where the las segment has j elements repeated r times.

                int CandidateLen = B[i-j*r].totalLen + j + rNDigits;
                if (CandidateLen < B[i].totalLen) {
                    B[i].lastSegment.elements = in + i - j*r;
                    B[i].lastSegment.nElements = j;
                    B[i].lastSegment.count = r;
                    B[i].prev = &(B[i-j*r]);
                    B[i].totalLen = CandidateLen;
                }

                r++;
                if (r == rNDigitsBound ) {
                    rNDigits++;
                    rNDigitsBound *= 10;
                }
            } while (   (i - j*r >= 0)
                     && (strncmp(in + i -j, in + i - j*r, j)==0));

        }
    }

    char *Res=sequence2string(&(B[N]));
    free(B);

    return Res;
}

int main(int argc, char** argv) {
    char *compressedDNA=dnaCompress(argv[1]);
    puts(compressedDNA);
    free(compressedDNA);
    return 0;
}

Забудь об Уконнене. Динамическое программирование это. С 3-мерной таблицей:

1. положение последовательности
2. размер подпоследовательности
3. количество сегментов

ТЕРМИНОЛОГИЯ: Например, при наличии a = "aaagctgctagag" координата позиции последовательности будет иметь значение от 1 до 13. В позиции последовательности 3 (буква 'g') при размере подпоследовательности 4 подпоследовательность будет "gctg". Понял? А что касается количества сегментов, то выражение a как "aaagctgctagag1" состоит из 1 сегмента (сама последовательность). Выражая его как «a3gct2ag2», он состоит из 3 сегментов. "aaagctgct1ag2" состоит из 2 сегментов. "a2a1ctg2ag2" будет состоять из 4 сегментов. Понял? Теперь, с этим, вы начинаете заполнять трехмерный массив 13 x 13 x 13, поэтому ваше время и сложность памяти, кажется, составляют около n ** 3 для этого. Вы уверены, что справитесь с последовательностями в миллион п.н.? Я думаю, что жадный подход был бы лучше, потому что большие последовательности ДНК вряд ли будут точно повторяться. И я бы посоветовал вам расширить свое задание до приблизительных совпадений, и вы можете опубликовать его прямо в журнале.

В любом случае, вы начнете заполнять таблицу сжатия подпоследовательности, начиная с некоторой позиции (размерность 1) с длиной, равной координате размерности 2, имеющей не более 3 сегментов размерности. Таким образом, вы сначала заполняете первую строку, представляющую сжатие подпоследовательностей длины 1, состоящих не более чем из 1 сегмента:

a        a        a        g        c        t        g        c        t        a        g        a        g
1(a1)    1(a1)    1(a1)    1(g1)    1(c1)    1(t1)    1(g1)    1(c1)    1(t1)    1(a1)    1(g1)    1(a1)    1(g1)

Число — это стоимость символа (всегда 1 для этих тривиальных 1-символьных последовательностей; число 1 не учитывается в стоимости символа), а в скобках у вас есть сжатие (также тривиальное для этого простого случая). Второй ряд будет по-прежнему прост:

2(a2)    2(a2)    2(ag1)   2(gc1)   2(ct1)   2(tg1)   2(gc1)   2(ct1)   2(ta1)   2(ag1)   2(ga1)    2(ag1)

Существует только 1 способ разложить 2-символьную последовательность на 2 подпоследовательности — 1 символ + 1 символ. Если они идентичны, результат похож на a + a = a2. Если они разные, например a + g, то, поскольку допустимы только последовательности из 1 сегмента, результат не может быть a1g1, а должен быть ag1. Третий ряд будет, наконец, более интересным:

2(a3)    2(aag1)  3(agc1)  3(gct1)  3(ctg1)  3(tgc1)  3(gct1)  3(cta1)  3(tag1)  3(aga1)  3(gag1)

Здесь вы всегда можете выбрать один из двух способов составления сжатой строки. Например, aag можно составить как aa + g, так и a + ag. Но опять же, у нас не может быть 2 сегментов, как в aa1g1 или a1ag1, поэтому мы должны довольствоваться aag1, если только оба компонента не состоят из одного и того же символа, как в aa + a => a3, со стоимостью символа 2. Мы можем продолжить на 4-м линия:

4(aaag1) 4(aagc1) 4(agct1) 4(gctg1) 4(ctgc1) 4(tgct1) 4(gcta1) 4(ctag1) 4(taga1) 3(ag2)

Здесь на первой позиции мы не можем использовать a3g1, потому что на этом слое разрешен только 1 сегмент. Но в последней позиции сжатие до стоимости символа 3 достигается за ag1 + ag1 = ag2. Таким образом, можно заполнить всю таблицу первого уровня вплоть до одной подпоследовательности из 13 символов, и каждая подпоследовательность будет иметь свою оптимальную стоимость символа и ее сжатие при ограничении первого уровня не более 1 связанного с ней сегмента. .

Затем вы переходите на 2-й уровень, где разрешено 2 сегмента... И снова снизу вверх вы определяете оптимальную стоимость и сжатие каждой координаты таблицы при ограничении количества сегментов данного уровня, сравнивая все возможные способы составьте подпоследовательность, используя уже вычисленные позиции, пока вы полностью не заполните таблицу и, таким образом, не вычислите глобальный оптимум. Есть некоторые детали, которые нужно решить, но извините, я не буду кодировать это для вас.

Попробовав некоторое время по-своему, я благодарен jbaylina за его прекрасный алгоритм и реализацию на языке C. Вот моя предпринятая версия алгоритма jbaylina в Haskell, а под ней дальнейшее развитие моей попытки алгоритма с линейным временем, который пытается сжимать сегменты, которые включают повторяющиеся шаблоны, один за другим:

import Data.Map (fromList, insert, size, (!))

compress s = (foldl f (fromList [(0,([],0)),(1,([s!!0],1))]) [1..n - 1]) ! n  
 where
  n = length s
  f b i = insert (size b) bestCandidate b where
    add (sequence, sLength) (sequence', sLength') = 
      (sequence ++ sequence', sLength + sLength')
    j' = [1..min 100 i]
    bestCandidate = foldr combCandidates (b!i `add` ([s!!i,'1'],2)) j'
    combCandidates j candidate' = 
      let nextCandidate' = comb 2 (b!(i - j + 1) 
                       `add` ((take j . drop (i - j + 1) $ s) ++ "1", j + 1))
      in if snd nextCandidate' <= snd candidate' 
            then nextCandidate' 
            else candidate' where
        comb r candidate
          | r > uBound                         = candidate
          | not (strcmp r True)                = candidate
          | snd nextCandidate <= snd candidate = comb (r + 1) nextCandidate
          | otherwise                          = comb (r + 1) candidate
         where 
           uBound = div (i + 1) j
           prev = b!(i - r * j + 1)
           nextCandidate = prev `add` 
             ((take j . drop (i - j + 1) $ s) ++ show r, j + length (show r))
           strcmp 1   _    = True
           strcmp num bool 
             | (take j . drop (i - num * j + 1) $ s) 
                == (take j . drop (i - (num - 1) * j + 1) $ s) = 
                  strcmp (num - 1) True
             | otherwise = False

Вывод:

*Main> compress "aaagctgctagag"
("a3gct2ag2",9)

*Main> compress "aaabbbaaabbbaaabbbaaabbb"
("aaabbb4",7)

Попытка линейного времени:

import Data.List (sortBy)

group' xxs sAccum (chr, count)
  | null xxs = if null chr 
                  then singles
                  else if count <= 2 
                          then reverse sAccum ++ multiples ++ "1"
                  else singles ++ if null chr then [] else chr ++ show count
  | [x] == chr = group' xs sAccum (chr,count + 1)
  | otherwise = if null chr 
                   then group' xs (sAccum) ([x],1) 
                   else if count <= 2 
                           then group' xs (multiples ++ sAccum) ([x],1)
                   else singles 
                        ++ chr ++ show count ++ group' xs [] ([x],1)
 where x:xs = xxs
       singles = reverse sAccum ++ (if null sAccum then [] else "1")
       multiples = concat (replicate count chr)

sequences ws strIndex maxSeqLen = repeated' where
  half = if null . drop (2 * maxSeqLen - 1) $ ws 
            then div (length ws) 2 else maxSeqLen
  repeated' = let (sequence,(sequenceStart, sequenceEnd'),notSinglesFlag) = repeated
              in (sequence,(sequenceStart, sequenceEnd'))
  repeated = foldr divide ([],(strIndex,strIndex),False) [1..half]
  equalChunksOf t a = takeWhile(==t) . map (take a) . iterate (drop a)
  divide chunkSize b@(sequence,(sequenceStart, sequenceEnd'),notSinglesFlag) = 
    let t = take (2*chunkSize) ws
        t' = take chunkSize t
    in if t' == drop chunkSize t
          then let ts = equalChunksOf t' chunkSize ws
                   lenTs = length ts
                   sequenceEnd = strIndex + lenTs * chunkSize
                   newEnd = if sequenceEnd > sequenceEnd' 
                            then sequenceEnd else sequenceEnd'
               in if chunkSize > 1 
                     then if length (group' (concat (replicate lenTs t')) [] ([],0)) > length (t' ++ show lenTs)
                             then (((strIndex,sequenceEnd,chunkSize,lenTs),t'):sequence, (sequenceStart,newEnd),True)
                             else b
                     else if notSinglesFlag
                             then b
                             else (((strIndex,sequenceEnd,chunkSize,lenTs),t'):sequence, (sequenceStart,newEnd),False)
          else b

addOne a b
  | null (fst b) = a
  | null (fst a) = b
  | otherwise = 
      let (((start,end,patLen,lenS),sequence):rest,(sStart,sEnd)) = a 
          (((start',end',patLen',lenS'),sequence'):rest',(sStart',sEnd')) = b
      in if sStart' < sEnd && sEnd < sEnd'
            then let c = ((start,end,patLen,lenS),sequence):rest
                     d = ((start',end',patLen',lenS'),sequence'):rest'
                 in (c ++ d, (sStart, sEnd'))
            else a

segment xs baseIndex maxSeqLen = segment' xs baseIndex baseIndex where
  segment' zzs@(z:zs) strIndex farthest
    | null zs                              = initial
    | strIndex >= farthest && strIndex > 0 = ([],(0,0))
    | otherwise                            = addOne initial next
   where
     next@(s',(start',end')) = segment' zs (strIndex + 1) farthest'
     farthest' | null s = farthest
               | otherwise = if start /= end && end > farthest then end else farthest
     initial@(s,(start,end)) = sequences zzs strIndex maxSeqLen

areExclusive ((a,b,_,_),_) ((a',b',_,_),_) = (a' >= b) || (b' <= a)

combs []     r = [r]
combs (x:xs) r
  | null r    = combs xs (x:r) ++ if null xs then [] else combs xs r
  | otherwise = if areExclusive (head r) x
                   then combs xs (x:r) ++ combs xs r
                        else if l' > lowerBound
                                then combs xs (x: reduced : drop 1 r) ++ combs xs r
                                else combs xs r
 where lowerBound = l + 2 * patLen
       ((l,u,patLen,lenS),s) = head r
       ((l',u',patLen',lenS'),s') = x
       reduce = takeWhile (>=l') . iterate (\x -> x - patLen) $ u
       lenReduced = length reduce
       reduced = ((l,u - lenReduced * patLen,patLen,lenS - lenReduced),s)

buildString origStr sequences = buildString' origStr sequences 0 (0,"",0)
   where
    buildString' origStr sequences index accum@(lenC,cStr,lenOrig)
      | null sequences = accum
      | l /= index     = 
          buildString' (drop l' origStr) sequences l (lenC + l' + 1, cStr ++ take l' origStr ++ "1", lenOrig + l')
      | otherwise      = 
          buildString' (drop u' origStr) rest u (lenC + length s', cStr ++ s', lenOrig + u')
     where
       l' = l - index
       u' = u - l  
       s' = s ++ show lenS       
       (((l,u,patLen,lenS),s):rest) = sequences

compress []         _         accum = reverse accum ++ (if null accum then [] else "1")
compress zzs@(z:zs) maxSeqLen accum
  | null (fst segment')                      = compress zs maxSeqLen (z:accum)
  | (start,end) == (0,2) && not (null accum) = compress zs maxSeqLen (z:accum)
  | otherwise                                =
      reverse accum ++ (if null accum || takeWhile' compressedStr 0 /= 0 then [] else "1")
      ++ compressedStr
      ++ compress (drop lengthOriginal zzs) maxSeqLen []
 where segment'@(s,(start,end)) = segment zzs 0 maxSeqLen
       combinations = combs (fst $ segment') []
       takeWhile' xxs count
         | null xxs                                             = 0
         | x == '1' && null (reads (take 1 xs)::[(Int,String)]) = count 
         | not (null (reads [x]::[(Int,String)]))               = 0
         | otherwise                                            = takeWhile' xs (count + 1) 
        where x:xs = xxs
       f (lenC,cStr,lenOrig) (lenC',cStr',lenOrig') = 
         let g = compare ((fromIntegral lenC + if not (null accum) && takeWhile' cStr 0 == 0 then 1 else 0) / fromIntegral lenOrig) 
                         ((fromIntegral lenC' + if not (null accum) && takeWhile' cStr' 0 == 0 then 1 else 0) / fromIntegral lenOrig')
         in if g == EQ 
               then compare (takeWhile' cStr' 0) (takeWhile' cStr 0)
               else g
       (lenCompressed,compressedStr,lengthOriginal) = 
         head $ sortBy f (map (buildString (take end zzs)) (map reverse combinations))

Вывод:

*Main> compress "aaaaaaaaabbbbbbbbbaaaaaaaaabbbbbbbbb" 100 []
"a9b9a9b9"

*Main> compress "aaabbbaaabbbaaabbbaaabbb" 100 []
"aaabbb4"