Как совместить шаблон последовательности dna

Ответ 1

Объяснение алгоритма:

• Имея последовательность 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), мы выберем наилучшую последовательность среди кандидатов:

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

  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) я = s (i + 1) i; только если s (1) = s (2) и s (2) = s (3) и... и s (i) = s (i + 1)

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

  В (I-1) с (я) с (г + 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 элемента:

  ...

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

  ...

  ...

• Последовательности кандидатов, где последний блок содержит n + 1 элементов:

  с (1) с (2) с (3)......... с (г + 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 строк кода.

EDIT Я немного изменил код, изменил имена переменных и добавил комментарий, чтобы быть более читаемым.

#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;
}

Ответ 2

Забыть Уконнена. Динамическое программирование. С трехмерной таблицей:

• положение последовательности
• размер подпоследовательности
• количество сегментов

ТЕРМИНОЛОГИЯ: Например, имея a = "aaagctgctagag", координата положения последовательности будет выполняться от 1 до 13. В позиции последовательности 3 (буква "g" ), имеющей размер 4 подпоследовательности, подпоследовательность будет "gctg". Понял? А что касается количества сегментов, то выражение "aaagctgctagag1" состоит из 1 сегмента (самой последовательности). Выражая это как "a3gct2ag2", состоит из 3 сегментов. "aaagctgct1ag2" состоит из 2 сегментов. "a2a1ctg2ag2" будет состоять из 4 сегментов. Понял? Теперь вы начинаете заполнять 3-мерный массив 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 char, номер 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)

Существует только один способ разложить 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 сегмента. И снова, снизу вверх, вы определяете оптимальную стоимость и сжатие каждой координаты таблицы под заданным ограничением счетчика уровня, сравнивая все возможные способы составления подпоследовательности с использованием уже вычисленных позиций, пока вы полностью не заполните таблицу и, следовательно, не вычислите глобальный оптимум. Есть некоторые детали, которые нужно решить, но, к сожалению, я не собираюсь писать это для вас.

Ответ 3

Попытка моего собственного пути на некоторое время, моя преданность 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 [email protected](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' [email protected](z:zs) strIndex farthest
    | null zs                              = initial
    | strIndex >= farthest && strIndex > 0 = ([],(0,0))
    | otherwise                            = addOne initial next
   where
     [email protected](s',(start',end')) = segment' zs (strIndex + 1) farthest'
     farthest' | null s = farthest
               | otherwise = if start /= end && end > farthest then end else farthest
     [email protected](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 [email protected](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 [email protected](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"