求教如何快速截取一个较长序列中的几段序列

feiyue0908

假设已知一个数组@allcdss=(2,8,12,16,18,25,28,35)；
已知一个文件，（只是举例，这个文件大约有上亿个字母），如下：
>chr5
AGAACAATATCTCTGATGAACATAGATGCAAAAATCCTCAACAACATACTAGCTAAATAATCCAACAACT
CATCAAAAGGATAATACACCATGATCAAGTGGGTTCATCCTAGAGATGCAGGGATGGTTTAACACATGCA
AGTCAATAAATATGATACATCACATAAACGGAATTAAAACAAAAACCACATGATCATCTCAATAGATGCA
GGAAAAGCATTTGATAAAATCTAGCAGCCCTTTATGATAAAACCCCTCAAGAAACTAGGTAAAGAAGGGG



想要的结果是根据数组内的值，截取中间的序列。2-8， 12-16， 18-25， 28-35,共四段，然后连接起来，
预计结果是GAACAATTCTGAGAACATAGGCAAAAAT



我用bioperl写过程序，把这个大文件读入就需要一个小时，要处理的数组和文件都多，这样耗费时间不大现实。
求教快速的方法，非常感谢

yinyuemi

回复 1# feiyue0908


    边读文件边截取的思路 应该可以完成你的要求

zhlong8

你内存多少，按你说的数据也就多一百多M怎么会读一个小时。

另外你这个可以按流的方式处理，一次读入一条数据截取然后就可以丢掉了

feiyue0908

回复 3# zhlong8


    上面的例子肯定很快就运行完，但完整的文件有1亿多个字母，用bioperl读入这个fasta文件就需要很长时间的。

cdtits

seek(fd, 2, 0);
read 8-2+1 bytes
seek(fd, 12, 0);
read 16-12+1 bytes;

其他依次类推，应该很快

zhlong8

feiyue0908 发表于 2012-07-27 09:15
回复 3# zhlong8

看来bioperl 做了很多额外的工作，但是这里明显影响性能了要么你就查文档这数据量应该在处理范围，要么就自己动手

kk861123

写了一个，你试试

1. #!/usr/bin/perl
   
2. use strict;
   
3. use warnings;
   
4. 
   
5. my @allcdss=(2,8,12,16,18,25,28,35,80,90,102,110,120,125);
   
6. my @allcuts;
   
7. {
   
8.     local $/ = \100;
   
9.     my ($min,$max) = (0,0);
   
10.     my $left = '';
    
11.     while (defined (my $line = <DATA>)) {
    
12.         $line =~ s/\n//g;
    
13.         $min  = $max - length($left);
    
14.         $max  = $max + length($line);
    
15.         $line = $left . $line;
    
16.         
    
17.         print "$min $max \n";
    
18.         my @allcdss_tmp  = grep { $_ > $min && $_ <= $max } @allcdss;
    
19.         if (@allcdss_tmp) {
    
20.             my $template = cut2fmt($min,@allcdss_tmp);
    
21.             print $template,"\n";
    
22.             my @fields   = unpack($template,$line);
    
23.             if ($allcdss_tmp[-1] != $max) {
    
24.                 $left = pop @fields;
    
25.             }
    
26.             push @allcuts,@fields;
    
27.         } else {
    
28.             $left = $line;
    
29.             next;
    
30.         }
    
31.         
    
32.         last if $allcdss_tmp[-1] == $allcdss[-1];
    
33.     }
    
34. }
    
35. print "@allcuts\n";
    
36. 
    
37. sub cut2fmt {
    
38.     my $minpos = shift;
    
39.     my(@positions) = @_;
    
40.     print "@positions\n";
    
41.     my $template   = ($positions[0] != $minpos) ? 'x' . ($positions[0] - $minpos - 1) : '';
    
42.     my $cut_string = 1;
    
43.    
    
44.     my $lastpos = shift @positions;
    
45.     foreach my $place (@positions) {
    
46.         if ($cut_string % 2 == 1) {
    
47.             my $dist   = $place - $lastpos;
    
48.             $template .= ($dist > 0) ? 'A' . ($dist + 1) : '';
    
49.         } else {
    
50.             my $dist   = $place - $lastpos;
    
51.             $template .= ($dist > 0) ? 'x' . ($dist - 1) : '';
    
52.         }
    
53.         $lastpos = $place;
    
54.         $cut_string++;
    
55.     }
    
56.       
    
57.     $template  .= 'A*';
    
58.     return $template;
    
59. }
    
60. 
    
61. __DATA__
    
62. AGAACAATATCTCTGATGAACATAGATGCAAAAATCCTCAACAACATACTAGCTAAATAATCCAACAACT
    
63. CATCAAAAGGATAATACACCATGATCAAGTGGGTTCATCCTAGAGATGCAGGGATGGTTTAACACATGCA
    
64. AGTCAATAAATATGATACATCACATAAACGGAATTAAAACAAAAACCACATGATCATCTCAATAGATGCA
    
65. GGAAAAGCATTTGATAAAATCTAGCAGCCCTTTATGATAAAACCCCTCAAGAAACTAGGTAAAGAAGGGG

复制代码