如何匹配dna序列模式

Question

我找不到解决这个问题的方法.

输入输出序列如下:

 **input1 :** aaagctgctagag 
 **output1 :** a3gct2ag2

 **input2 :** aaaaaaagctaagctaag 
 **output2 :** a6agcta2ag

输入序列可以是10 ^ 6个字符,并且将考虑最大的连续模式.

例如,对于input2"agctaagcta"输出将不是"agcta2gcta",但它将是"agcta2".

任何帮助赞赏.

Answer 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),我们将在候选人中选择最佳序列:

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

Answer 2

忘记乌科南吧。就是动态规划。使用 3 维表：

1. 序列位置
2. 子序列大小
3. 段数

术语：例如，a = "aaagctgctagag"序列位置坐标将从 1 到 13。在序列位置 3（字母“g”）处，子序列大小为 4，子序列将为“gctg”。明白了吗？至于分段数，则将 a 表示为“aaagctgctag1”，由 1 个分段（序列本身）组成。将其表示为“a3gct2ag2”，由 3 段组成。“aaagctgct1ag2”由 2 段组成。“a2a1ctg2ag2”将由 4 个段组成。明白了吗？现在，您开始填充一个 3 维数组 13 x 13 x 13，因此您的时间和内存复杂度似乎n ** 3与此相关。您确定可以处理百万bp的序列吗？我认为贪婪的方法会更好，因为大的 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)

只有 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(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 个段...并且，从下到上，您可以通过比较所有可能的方法来确定给定级别的段计数约束下每个表坐标的最佳成本和压缩使用已经计算出的位置组成子序列，直到完全填满表格并计算出全局最优值。有一些细节需要解决，但抱歉，我不会为您编写代码。