[求助]在5万个长3k的字符串中快速定位定长字符串

galaxy001

UID: 40457
帖子: 169
积分: 388
在线时间: 1 天 15 小时

1^# galaxy001 发表于 2007-09-20 17:12

[求助]在5万个长3k的字符串中快速定位定长字符串

我需要做批字符串定位，目标数据文件格式为：

[Copy to clipboard] [ - ]

CODE:

>字符串名
字符串内容
>下一个字符串名
字符串内容

以下为实际例子，字符串长度有删节：

[Copy to clipboard] [ - ]

CODE:

>LOC_Os01g01010.1|12001.m06748|cDNA TBC domain containing protein, expressed
CTCCCCGCCCACCAGGCCGCCAGCCACAACCAGCTGACTAGGTCGCATGCATCATCAGAT
TTCAATCTCCCTTCGTTCCCTGTCCCTAATCCAATACCAATAGGGAGCAATCAGCTGCTC
>LOC_Os01g01010.2|12001.m42814|cDNA TBC domain containing protein, expressed
ACGACGCAGCTATGGCCTCCCCGCCCACCAGGCCGCCAGCCACAACCAGCTGACTAGGTA
GGCTTCCTAGGTCGCATGCATCATCAGATTTCAATCTCCCTTCGTTCCCTGTCCCTAATC
>LOC_Os08g20470.1|12008.m06113|cDNA retrotransposon protein, putative, unclassified
ATGGTTGATGATTGTTTTAAGTACAATAGAGGATGTGAAGCTTGTCAACAATTCGGCAAT

单一文件，大小为120M。
每个字符串长3k左右，构成字符串的字符只有A,T,C,G四种。

待检索字符串恒为25字节，检索结果是待检索字符串在哪个目标数据字符串的第几字节。可能有多个匹配结果。
如GACGCAGCTATGGCCTCCCCGCCCA在>LOC_Os01g01010.2|12001.m42814|cDNA TBC domain containing protein, expressed的第3位。

我想了一个方法，但没能力写代码：
每目标字符串连续每25字节(1..25,2..26,3..27,...)做crc8的hash，对hash做索引（好像一般用btree?），查询时把待检索字符串算crc8，对匹配项逐个作普通字符串比较(eq)。
至于hash长度，我也不知道多长合适。25字节的信息熵是50bit，好像crc4就可以了？长了费内存，短了每个hash对应的字符串太多影响时间效率。

能写个较完整的代码最好，只有搜索部分我也会满足的。

Lonki

UID: 21498
帖子: 1
积分: 2
在线时间: 10 分钟

2^# Lonki 发表于 2007-09-21 00:15

"每个字符串长3k左右"

从例子的格式看, 这个3K的ACGT字符串似乎由于多行组成....

galaxy001

UID: 40457
帖子: 169
积分: 388
在线时间: 1 天 15 小时

3^# galaxy001 发表于 2007-09-21 01:00

所以用>或FEOF作为上一个字符串的结束。

perljoker

UID: 30669
帖子: 88
积分: 202
在线时间: 10 小时

4^# perljoker 发表于 2007-09-21 10:54

QUOTE:

原帖由 galaxy001 于 2007-9-21 01:00 发表
所以用>或FEOF作为上一个字符串的结束。

Lonki的意思应该是指你的fasta序列的dna序列是多行文本
如果是单行就好的多了$_=~/$string/ig
位置可以用pos()得到，可以参考perldoc
以前我用顺序读取待测序列的方法，相当低效
正则表达式的合理使用很高效
关键是你如何处理多行文本

Lonki

UID: 21498
帖子: 1
积分: 2
在线时间: 10 分钟

5^# Lonki 发表于 2007-09-21 13:00

QUOTE:

原帖由 perljoker 于 2007-9-21 10:54 发表

Lonki的意思应该是指你的fasta序列的dna序列是多行文本
如果是单行就好的多了$_=~/$string/ig
位置可以用pos()得到，可以参考perldoc
以前我用顺序读取待测序列的方法，相当低效
正则表达式的合理使用很高 ...

对, 正是这意思. 不过对于100来M的文件, 还是用Tie::File吧, 用LZ的例子作为输入文件.
此处用index比正则更高效... 正则用习惯了 :>

tie my @lines, 'Tie::File', $ARGV[0] or die "Where the hell is $ARGV[0]?";

    my ($pattern, $name, $str) = ('GA', '', '');

    for (@lines) {
        if (substr($_, 0, 1) eq '>') {
            while ($str =~ /$pattern/g) {
                print "Found pattern in pos ",  $-[0] + 1, " of $name\n";
            }
            $name = $_;
            $str  = '';
        } else {
            $str .= $_;
        }
    }

    untie @lines;
}

------------------- Running Result ------------------------
Found pattern in pos 36 of >LOC_Os01g01010.1|12001.m06748|cDNA TBC domain containing protein, expressed
Found pattern in pos 58 of >LOC_Os01g01010.1|12001.m06748|cDNA TBC domain containing protein, expressed
Found pattern in pos 105 of >LOC_Os01g01010.1|12001.m06748|cDNA TBC domain containing protein, expressed
Found pattern in pos 3 of >LOC_Os01g01010.2|12001.m42814|cDNA TBC domain containing protein, expressed
Found pattern in pos 52 of >LOC_Os01g01010.2|12001.m42814|cDNA TBC domain containing protein, expressed
Found pattern in pos 87 of >LOC_Os01g01010.2|12001.m42814|cDNA TBC domain containing protein, expressed

上面的代码作废....
Updated: 09/26/2007 11:23 AM
11楼 Perl 代码
16楼 C 代码

galaxy001

UID: 40457
帖子: 169
积分: 388
在线时间: 1 天 15 小时

6^# galaxy001 发表于 2007-09-23 13:54

搜一个要307s，刨掉tie文件的时间，也不行啊。

$-[0]是哪的用法？我翻大骆驼找到file::handle去了。

perljoker

UID: 30669
帖子: 88
积分: 202
在线时间: 10 小时

7^# perljoker 发表于 2007-09-23 16:42

QUOTE:

原帖由 galaxy001 于 2007-9-23 13:54 发表
搜一个要307s，刨掉tie文件的时间，也不行啊。

$-[0]是哪的用法？我翻大骆驼找到file::handle去了。

你一个文件有多大？

Lonki

UID: 21498
帖子: 1
积分: 2
在线时间: 10 分钟

8^# Lonki 发表于 2007-09-23 18:35

QUOTE:

原帖由 galaxy001 于 2007-9-23 13:54 发表
搜一个要307s，刨掉tie文件的时间，也不行啊。

$-[0]是哪的用法？我翻大骆驼找到file::handle去了。

1. 兄弟, 好歹也是100M的文件..... (能告知Tie的耗时吗?)
"搜一个"是指搜一个模式在整个文件?

2. 看本版精华之特殊变量
$-[0]表示上一个匹配到的字符串在原字符串中的起始位置.

galaxy001

UID: 40457
帖子: 169
积分: 388
在线时间: 1 天 15 小时

9^# galaxy001 发表于 2007-09-23 21:20

搜一个当然是用搜一个模式在整个文件，因为有可能有多个匹配，本来就没考虑中途跳出的。
Tie文件的时间好像没法单独测：

[Copy to clipboard] [ - ]

CODE:

#!/usr/bin/perl -w
use strict;
use warnings;
use Time::HiRes qw ( gettimeofday tv_interval );
use Tie::File;
my $start_time = [gettimeofday];
tie my @lines, 'Tie::File', $ARGV[0] or die "Where the hell is $ARGV[0]?";
my $stop_time = [gettimeofday];
print STDERR "\n Time1 Elapsed:\t",tv_interval( $start_time, $stop_time )," second(s)\n";
my $start_time2 = [gettimeofday];
my ($pattern, $name, $str) = ('TTCAATCTCCCTTCGTTCCCTGTCC', '', '');
my $pattern2='AGGGCACTAGACAGGCTGCAGCAAC';
for (@lines) {
      if (substr($_, 0, 1) eq '>') {
         while ($str =~ /$pattern/g) {
            # print "Found pattern in pos ",  $-[0] + 1, " of $name\n";
            print $name,"\t",$-[0] + 1,"\n";
         }
         $name = $_;
         $str  = '';
      } else {
         $str .= $_;
      }
}
my $stop_time2 = [gettimeofday];
print STDERR "\n Time2 Elapsed:\t",tv_interval( $start_time2, $stop_time2 )," second(s)\n";
my $start_time3 = [gettimeofday];
      for (@lines) {
      if (substr($_, 0, 1) eq '>') {
         while ($str =~ /$pattern2/g) {
            # print "Found pattern in pos ",  $-[0] + 1, " of $name\n";
            print $name,"\t",$-[0] + 1,"\n";
         }
         $name = $_;
         $str  = '';
      } else {
         $str .= $_;
      }
}
untie @lines;
my $stop_time3 = [gettimeofday];
print STDERR "\n Time3 Elapsed:\t",tv_interval( $start_time2, $stop_time2 )," second(s)\n";

结果是：

Time1 Elapsed: 0.000157 second(s)
>LOC_Os01g01010.1|12001.m06748|cDNA TBC domain containing protein, expressed 61
>LOC_Os01g01010.2|12001.m42814|cDNA TBC domain containing protein, expressed 90

Time2 Elapsed: 308.979046 second(s)
>LOC_Os09g36990.1|12009.m06648|cDNA expressed protein 412

Time3 Elapsed: 308.979046 second(s)

galaxy001

UID: 40457
帖子: 169
积分: 388
在线时间: 1 天 15 小时

10^# galaxy001 发表于 2007-09-23 23:16

for (@lines) {

$name = $_;

}
光这个循环就309s
开来瓶颈在文件IO，……