寻找连续重复序列的算法

Question

我正在寻找一种在基因组序列中发现短串联重复的算法.

基本上,给定一个非常长的字符串,只能包含4个字符'ATCG',我需要找到长度在2-5个字符之间的短重复,它们彼此相邻.

例如:TACATGAGATCATGATGATGATGATGGAGCTGTGAGATC将 重复ATGATGATG或ATG 3次

该算法需要扩展到一百万个字符的字符串,所以我试图尽可能接近线性运行时.

我的当前算法:由于重复可以是2-5个字符长,我通过字符检查字符串字符,看看是否第N个字符相同的N + X个字符,其中X为2至5.对于每个计数器X计算顺序匹配并在不匹配时重置,我们知道当X =计数器时是否重复.然后可以手动检查随后的重复.

Answer 1

您正在查看给您提供的每个字符O(n)，因为您比较每个字符接下来（最多）五个字符，这给您一个常数c：

var data    = get_input();
var compare = { `A`, `T`, `G`, `A`, `T` }         // or whatever
var MAX_LOOKAHEAD = compare.length
var n
var c

for(n = data_array.length; n < size; i++) {       // Has runtime O(n)

  for(c = 0; c < MAX_LOOKAHEAD; c++) {            // Maximum O(c)

    if( compare[c] != data[i+c] ) {
      break;
    } else {
      report( "found match at position " + i )
    }

  }
}

很容易看出这个运行O(n*c)时间。由于c非常小，因此可以忽略 - 而且我认为无法摆脱该常量 - 这导致总运行时间为O(n).

好消息：

您可以通过并行化来加快速度。例如，您可以将其分成一定的k时间间隔，并通过为多个线程提供适当的开始和结束索引来让它们为您完成这项工作。这可以给你带来线性加速。

如果您这样做，请确保将交叉点视为特殊情况，因为如果您的间隔将一场比赛分成两部分，您可能会错过一场比赛。

例如n = 50000：

4 个线程的分区：(n/10000) - 1 = 4. 第 5 个线程不会有很多事情要做，因为它只处理交叉点，这就是为什么我们不需要考虑它的（在我们的例子中很小）开销。

1                 10000               20000               40000               50000
|-------------------|-------------------|-------------------|-------------------|
| <-   thread 1  -> | <-   thread 2  -> | <-   thread 3  -> | <-   thread 4  -> |
                  |---|               |---|               |---|              
                    |___________________|___________________|
                                        |
                                     thread 5

它可能是这样的：

var data;
var compare = { `A`, `T`, `G`, `A`, `T` };
var MAX_LOOKAHEAD = compare.length;

thread_function(args[]) {

    var from = args[0];
    var to   = args[1];

    for(n = from ; n < to ; i++) {

      for(c = 0; c < MAX_LOOKAHEAD; c++) {
        if( compare[c] != data[i+c] ) {
          break;
        } else {
          report( "found match at position " + i )
        }
      }
    }
}

main() {
    var data_size     = 50000;
    var thread_count  = 4;
    var interval_size = data_size / ( thread_count + 1) ;

    var tid[]

    // This loop starts the threads for us:

    for( var i = 0; i < thread_count; i++ ) {
        var args = { interval_size * i, (interval_size * i) + interval_size };

        tid.add( create_thread( thread_function, args ) );
    }

    // And this handles the intersections:

    for( var i = 1; i < thread_count - 1; i++ ) {
        var args = { interval_size * i, (interval_size * i) + interval_size };

        from = (interval_size * i) - compare.length + 1;
        to   = (interval_size * i) + compare.length - 1;

        for(j = from; j < to ; j++) {

            for(k = 0; k < MAX_LOOKAHEAD; k++) {
                if( compare[k] != data[j+k] ) {
                    break;
                } else {
                    report( "found match at position " + j )
                }
            }
        }
    }

    wait_for_multiple_threads( tid );
}