如何匹配 dna 序列模式
Posted
技术标签:
【中文标题】如何匹配 dna 序列模式【英文标题】：how to match dna sequence pattern 【发布时间】：2013-05-28 02:11:24 【问题描述】：我很难找到解决这个问题的方法。
输入输出顺序如下：
**input1 :** aaagctgctagag
**output1 :** a3gct2ag2
**input2 :** aaaaaaagctaagctaag
**output2 :** a6agcta2ag
输入 nsequence 可以是 10^6 个字符，并且将考虑最大的连续模式。
例如，对于 input2，“agctaagcta”输出将不是“agcta2gcta”，而是“agcta2”。
任何帮助表示赞赏。
【问题讨论】：
输入aabbaabb必须提供什么输出？两种可能的变体：a2b2a2b2 和 aabb2。
输出应该是“aabb2”
那么aaaaaaaaabbbbbbbbbaaaaaaaaabbbbbbbbb：a9b9a9b9 或aaaaaaaaabbbbbbbbb2 呢？前者更短;-)
字符数及其计数应该是最少的..例如 a9b9a9b9 需要 8 个字母数字，但 aaaaaaaaaabbbbbbbbb2 需要 19 个字母数字
你会如何编码这个：aaagctgctxyzagag?
【参考方案1】：
算法说明：
有一个带有符号 s(1)、s(2)、...、s(N) 的序列 S。 令 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) and s(3)=s(4) and … and s(i) = s(i+1) B(0)s(i+1)i = s(i+1)i ;仅当 s(1)=s(2) and s(2)=s(3) and … and 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 个元素的候选序列：
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;
【讨论】：
+1，似乎工作得很好，C 也是一种很好的语言来实现这一点，虽然有点难以阅读。例如。ss 变量是什么。这个解决方案是最优的吗？给他一个这样的解决方案对 OP 来说有点太容易了。他应该在他的论文中给你荣誉:-)))
看到你现在的名声和应得的名声之间的差异，一定要为此赏金。
我保证会尝试多解释一点算法。无论如何，s 是 r 的位数，ss 是将 r 的密码数加一的下一个数字。
算法应该给出最优值。如果您想限制可以重复的最大序列长度，如您在其中一个 cmets 中所说，您可以限制第二个 for 循环运行最大 N（其中 N 是最大序列）。所以该行将是： for (j=i-i, (j>0)&&(i-j
@jbaylina：所以它毕竟是动态编程。到目前为止，我一生中写的 C 很少，无法理解您的代码。这种问题和答案是什么让如此有趣。在您解释了魔术之后，这听起来几乎像天真的算法:-)。但是，我想知道，这可以在较短的 O 中完成吗？比如n log n或者linear，用后缀树什么的？【参考方案2】：
忘记乌克南吧。动态规划就是这样。带3维表：
-
序列位置
子序列大小
段数
术语：例如，有a = "aaagctgctagag"，序列位置坐标将从1 到13。在序列位置3（字母'g'），子序列大小为4，子序列将是“gctg”。明白了吗？至于段数，则将a表示为“aaagctgctagag1”由1段（序列本身）组成。将其表示为“a3gct2ag2”由 3 个部分组成。 “aaagctgct1ag2”由 2 个段组成。 “a2a1ctg2ag2”将由 4 个段组成。明白了吗？现在，有了这个，您开始填充一个 13 x 13 x 13 的 3 维数组，因此您的时间和内存复杂度似乎在 n ** 3 左右。你确定你可以处理百万碱基序列吗？我认为贪婪的方法会更好，因为大的 DNA 序列不太可能完全重复。而且，我建议您将作业范围扩大到近似匹配，并且可以直接在期刊上发表。
无论如何，您将开始填充从某个位置（维度 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)
只有一种方法可以将 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)
在这里，您始终可以在 2 种组合压缩字符串的方式中进行选择。例如，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，因为在这一层只允许有一个段。但在最后一个位置，ag1 + ag1 = ag2 实现了对字符成本 3 的压缩。这样，一个人可以将整个一级表一直填充到 13 个字符的单个子序列，并且每个子序列将在与其关联的最多 1 个段的一级约束下具有其最佳字符成本和压缩.
然后进入第二级，允许 2 个段...再次，从下往上，通过比较所有使用已经计算的位置组成子序列的可能方法，直到您完全填满表格并因此计算全局最优值。有一些细节需要解决，但抱歉，我不会为你编写代码。
【讨论】：
这大约是 O(n^4)，因为每个元素的计算都包含一个内部循环。 @EgorSkriptunoff：这也许可以避免。也许你甚至可以在 n^2 中实现它，如果你不娇气，你只需检查子序列达到一定长度（比如说 4000 bp），然后保存指数。但是我已经花了太多时间思考这个问题。这就是为什么我提到有一些细节需要解决。这个问题实际上具有可疑的线性结构，我怀疑这个 O(n^3) （或者可能是 n^4，如你所说）存在更好的解决方案。问题只是找到一个足够好的。指数时间不行。 @EgorSkriptunoff：但是是的，我必须把它给你，它似乎是n ** 4，因为必须比较多个组合，或者甚至更高的指数，最多 5，我'我真的懒得数了。【参考方案3】：
在尝试了一段时间后，我对 jbaylina 的漂亮算法和 C 实现表示敬意。这是我在 Haskell 中尝试的 jbaylina 算法版本，下面是我对线性时间算法的尝试的进一步开发，该算法尝试以一个接一个的方式压缩包含重复模式的片段：
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"
【讨论】：
“aaabbbaaabbbaaabbbaaabbb”怎么样，正确答案是“aaabbb4”？ @BorisStitnicky 感谢您注意到这一点。我将divided = foldr divide [] [div lenOne 2, div lenOne 2 - 1..2] 更改为divided = foldr divide [] [2..div lenOne 2]，现在它显示aaabbb4。但这种变化会导致其他不准确吗？
@BorisStitnicky ..它比我想象的要复杂..我又试了一下
你不必告诉我。我贪婪地尝试用 Ruby 编写它，但是没有动态编程的计算复杂度实在是太糟糕了。由于问题不是教科书作业，所以缓慢的解决方案并不重要，所以我放弃了。
@BorisStitnicky 我做了一点修改。现在，它会在大约一分钟内将随机 10^6 字符串的压缩版本输出到文件中。它应该更快吗？以上是关于如何匹配 dna 序列模式的主要内容，如果未能解决你的问题，请参考以下文章
POJ2778 DNA Sequence Trie+矩阵乘法
python DNA基序匹配操作 - 查找与TATAA序列重叠的峰