Как соответствовать шаблону последовательности днк
У меня проблемы с поиском подхода для решения этой проблемы.
Последовательности ввода-вывода следующие:
**input1 :** aaagctgctagag
**output1 :** a3gct2ag2
**input2 :** aaaaaaagctaagctaag
**output2 :** a6agcta2ag
Входная последовательность может содержать 10^6 символов, и будут рассматриваться самые большие непрерывные шаблоны.
Например, для input2 "agctaagcta" вывод не будет "agcta2gcta", но будет "agcta2".
Любая помощь приветствуется.
3 ответа
Объяснение алгоритма:
- Имея последовательность 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)i = s(i+1) я; только если s(1)=s(2) и s(2)=s(3) и… и s (i) = s (i + 1)
Последовательности кандидатов, где последний блок имеет 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 элемента:
...
Последовательности кандидатов, где последний блок имеет 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 строк кода.
РЕДАКТИРОВАТЬ Я немного реструктурировал код, изменил имена переменных и добавил некоторые комментарии, чтобы быть более читабельным.
#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;
}
Попробовав свой собственный путь некоторое время, я благодарю 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"
Забудь Ukonnen. Динамическое программирование это так. С 3-х мерным столом:
- позиция последовательности
- размер подпоследовательности
- количество сегментов
ТЕРМИНОЛОГИЯ: Например, имея 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 символа; число 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 сегмента... И снова снизу вверх вы определяете оптимальную стоимость и сжатие каждой координаты таблицы в соответствии с ограничением количества сегментов данного уровня, сравнивая все возможные способы составьте подпоследовательность, используя уже вычисленные позиции, пока вы не заполните таблицу полностью и, таким образом, не вычислите глобальный оптимум. Есть некоторые детали, которые нужно решить, но извините, я не собираюсь кодировать это для вас.