perl的内存占用，求教&探讨

chenhao392

今天有个生物信息的任务，完成了，但是很不爽的搞死了电脑两次。于是全程播报我的思路和步骤，来论坛求教了，求拍砖：

1. 数据描述

有基因组数据如下，大写的是正常序列，小写的是repeat区域：

1. >chr1
   
2. GAATTCCAAAGCCAAAGATTGCATCAGTTCTGCTGCTATTTCCTCCTATCATTCTTTCTG
   
3. ATGTTGAAAATGATATTAAGCCTAGGATTCGTGAATGGGAGAAGGTATTTTTGTTCATGG
   
4. TAGTCATTGGAACCTGCTAGATTGTACACTTGACAATAACATATATTAATATTAGTGACC
   
5. CCATTTTTAAATTTCCTAGGCTGGCATTGAACAAGACTATGTTAGTAGGATGTTGTTGAA
   
6. GTATCCATGGATTCTTTCAACGAGTGTGATAGAGAACTACAGTCAAATGCTGTTGTTTTT
   
7. CAACCAAAAAAGGGTAAGTAAAAAAGAATACTTACTATGCTGTGCCTCAAGTTCATGTTA
   
8. 
   
9. TATTCAAATGCCGCAGCTCTGAtaaccactcttttctggaccaataaatggctgcttggt
   
10. tcctctataggagatgtgtcgctggctgttcttgctatggtacccgcaaaaataattaaa
    
11. 
    
12. cgtcttcgaAGGTTTGCAGCTGAGTATGTGGAATGTTCTGCATGCTCTCGAAGATGAACA
    
13. GCATCTCTGGTGTCTCGCCGGAGTGCAAGGACTGCAGAGTTTAGGCCTGGGCTCCCTGAG
    
14. CGGTTAGGTGGCATTTCCTGGTCTGGTCTCACCATTTCTATTTGCTGTAAGAGTTTTGTT
    
15. >chr2
    
16. .............
    
17. .............
    

复制代码

2. 任务要求
按照小写的序列区域，将一个长的染色体分割成较短的片段，小写区域丢弃，大写片段短于500字符就丢弃。

3. 分析每一行的可能格式

1.
>chr1

2.
大写片段-----小写片段-----大写片段
小写----大写---小写（由于每行60字符，这里的大写片段显然不够500字符）
小写------大写
大写------小写

3.
一直小写
一直大写

4. 写代码如下：

1. #!/usr/bin/perl -w
   
2. #
   
3. #Author **
   
4. #This script is uesd for matching the repeat masked region, and then split
   
5. #by them.
   
6. 
   
7. use POSIX;
   
8. 
   
9. my $repeat=$ARGV[0];#file repeat.fa
   
10. my $out=$ARGV[1];#splited file
    
11. my $min_seq_length=500;
    
12. my $out_line_length=50;
    
13. 
    
14. #repeat file load
    
15. open REPEAT,"<$repeat" or die "hi:$!";
    
16. my %chrom;
    
17. my $base;
    
18. while(<REPEAT>){
    
19.         chomp;
    
20.         if($_ =~ /^\>(.*?)$/){
    
21.                 $base=$1;
    
22.                 $chrom{$base}="";
    
23.         }
    
24.         else{
    
25.                 $chrom{$base}.=$_;
    
26.         }
    
27. }
    
28. close REPEAT;
    
29. 
    
30. #print file splited by repeats
    
31. open OUT,">$out" or die"Oh:$!";
    
32. foreach $base (sort keys %chrom){
    
33. 
    
34.         my @seq=split(//,$chrom{$base});
    
35.         my $count=0;
    
36.         my $start;
    
37.         my $stop;
    
38.         my $seq="";
    
39.         my $switch=0;#1 == ATCG, 0 == N
    
40. 
    
41.         foreach $a(@seq){
    
42.                 if($a !~ /[atcg]/ && $switch == 0){  #new seq after repeat
    
43.                         $switch=1;
    
44.                         $start=$count+1;
    
45.                         $seq.=$a;
    
46.                 }
    
47.                 elsif($a !~/[atcg]/ && $switch == 1){
    
48.                         $seq.=$a;
    
49.                 }
    
50.                 elsif($a =~ /[atcg]/ && $switch == 1){
    
51.                         $switch=0;
    
52.                         $stop=$count;
    
53. 
    
54.                         if(  ($stop - $start)   >=    $min_seq_length){
    
55.                                 print OUT "\>$base\_$start\_$stop\n";
    
56.                                 seq_print($seq,$out_line_length);
    
57.                         }
    
58. 
    
59.                         $seq="";
    
60.                 }
    
61. 
    
62.                 $count++;
    
63.         }
    
64.         #print the last seq
    
65.         if($switch == 1){$stop=$count-1;}
    
66.         if(     ($stop - $start)   >=    $min_seq_length){
    
67.                 print OUT "\>$base\_$start\_$stop\n";
    
68.                 seq_print($seq,$out_line_length);
    
69.         }
    
70. 
    
71. print "work on $base\n";
    
72. }
    
73. close OUT;
    
74. 
    
75. 
    
76. 
    
77. sub seq_print{
    
78.         my ($seq,$length)=@_;
    
79.         my $step=floor(length($seq)/$length);
    
80.         for(my $i=1;$i<=$step;$i++){
    
81.                 my $out=substr($seq,($i-1)*$length,$length);
    
82.                 print OUT "$out\n";
    
83.         }
    
84.         my $out=substr($seq,$step*$length);
    
85.         print OUT "$out\n";
    
86. }

复制代码

5. 测试小量数据通过
......


6. 大数据量，死机
于是在应用在了基因组上，数据量2G多

瞬间占用40G以上内存，服务器死机

7. 修改代码
想了想，我的@seq可能太过于恐怖了，于是修改代码如下：

1. #!/usr/bin/perl -w
   
2. #
   
3. #Author **
   
4. #This script is uesd for matching the repeat masked region NNNNN, and then split
   
5. #by them.
   
6. 
   
7. use POSIX;
   
8. 
   
9. my $repeat=$ARGV[0];#file repeat.fa
   
10. my $out=$ARGV[1];#splited file
    
11. my $min_seq_length=500;
    
12. my $out_line_length=50;
    
13. 
    
14. #repeat file load
    
15. open REPEAT,"<$repeat" or die "hi:$!";
    
16. my %chrom;
    
17. my $base;
    
18. my $count=0;
    
19. while(<REPEAT>){
    
20.         chomp;
    
21.         if($_ =~ /^\>(.*?)$/){
    
22.                 $base=$1;
    
23.                 $count=0;
    
24.         }
    
25.         else{
    
26.                 $chrom{$base}[$count]=$_;
    
27.                 $count++;
    
28.         }
    
29. }
    
30. close REPEAT;
    
31. 
    
32. 
    
33. open OUT,">$out" or die"Oh:$!";
    
34. foreach $base (sort keys %chrom){
    
35. 
    
36.         my $array_hash=$chrom{$base};
    
37.         my $count=0;
    
38.         my $start;
    
39.         my $stop;
    
40.         my $seq="";
    
41.         my $switch=0;#1 == ATCG, 0 == N
    
42. 
    
43.         foreach $name(@$array_hash){
    
44.           my @line=split(//,$name);
    
45.           foreach $a(@line){
    
46.                 if($a !~ /[atcg]/ && $switch == 0){  # first base of seq after repeat
    
47.                         $switch=1;
    
48.                         $start=$count+1;
    
49.                         $seq.=$a;
    
50.                 }
    
51.                 elsif($a !~ /[atcg]/ && $switch == 1){ #continue seq
    
52.                         $seq.=$a;
    
53.                 }
    
54.                 elsif($a =~ /[atcg]/ && $switch == 1){ #come across repeat
    
55.                         $switch=0;
    
56.                         $stop=$count;
    
57. 
    
58.                         if(  ($stop - $start)   >=    $min_seq_length){
    
59.                                 print OUT "\>$base\_$start\_$stop\n";
    
60.                                 seq_print($seq,$out_line_length);
    
61.                         }
    
62. 
    
63.                         $seq="";
    
64.                 }
    
65. 
    
66.                 $count++;
    
67.           }
    
68.         undef(@line);
    
69. 
    
70.         }
    
71.         #print the last seq
    
72.         if($switch == 1){$stop=$count-1;}
    
73.         if(     ($stop - $start)   >=    $min_seq_length){
    
74.                 print OUT "\>$base\_$start\_$stop\n";
    
75.                 seq_print($seq,$out_line_length);
    
76.         }
    
77. 
    
78. print "work on $base\n";
    
79. }
    
80. close OUT;
    
81. 
    
82. 
    
83. 
    
84. sub seq_print{
    
85.         my ($seq,$length)=@_;
    
86.         my $step=floor(length($seq)/$length);
    
87.         for(my $i=1;$i<=$step;$i++){
    
88.                 my $out=substr($seq,($i-1)*$length,$length);
    
89.                 print OUT "$out\n";
    
90.         }
    
91.         my $out=substr($seq,$step*$length);
    
92.         print OUT "$out\n";
    
93. }

复制代码

8. 求教
这次占用了大概4.3G的内存，程序顺利跑完，但是我依然不明白哪里需要4G内存了。
在我看来，这两行：
$seq="";
undef(@line);
已经释放了内存，应该只是%chrom存储的2G序列而已了。

zhlong8

2G 的原始数据才占 4G 的内存，已经很省了啊？

chenhao392

回复 2# zhlong8


    这个....这样啊....
    那另外2G内存被什么吃掉了？

    或者说，
    为什么4G就很省了？

zhlong8

将一个长的染色体分割成较短的片段，小写区域丢弃，大写片段短于500字符就丢弃。

这个用 RE 就能做了啊？你为什么一个一个字符处理？

@seq 那里 split 成单个字符每1K个占用内存是 48K，如果你一段基因超长的话……

Perl 的字符串可变，所以可能分配的内存比需要的多，而且是以2^n 的速度扩大的，所以平均长度越长浪费的也就越多

数据结构不也要占用内存，建一个空空的数组或 hash 难道就不需要空间了？

chenhao392

回复 4# zhlong8


谢谢说明！


将一个长的染色体分割成较短的片段，小写区域丢弃，大写片段短于500字符就丢弃。
这个用 RE 就能做了啊？你为什么一个一个字符处理？

我只能说是，想到什么就写什么了，另，RE是什么？


@seq 那里 split 成单个字符每1K个占用内存是 48K，如果你一段基因超长的话……

是的，恐怖的长，2G只有10个基因组，平均一个200MB. 看来这就是我第一次程序挂掉的原因


Perl 的字符串可变，所以可能分配的内存比需要的多，而且是以2^n 的速度扩大的，所以平均长度越长浪费的也就越多
数据结构不也要占用内存，建一个空空的数组或 hash 难道就不需要空间了？

$seq的长度经常会有几千上万，一直在变，那这就是原因了？
我是不是可以理解为：

$a="A";
$b="B";
$c="CCCCC";


$a="AA";
$b="BB";
$c="CC";

前者比后者 消耗的内存多？

zhlong8

对于很长的数据按行存储和处理比较高效。

说用 RE 指的是对于每一行的数据可以用正则表达式来判断分割。

对于比较短的字符串不需要考虑这么多，一个变量保存字符串的长度是曾经最长的长度，即使赋值了个比较短的字符串，相当于剩余的空间就是缓冲区了。用 undef 可以释放

chenhao392

回复 6# zhlong8


谢谢，变量占用内存大概明白了...


    RE指，类似这样的处理么？
     if (   $line=~/([ATCG]*)([atcg]*)([ATCG]*)/ ) {
                $1;
                $2;
                $3;
   }

其实一开始我是有想这么写的，但是还需要计算大写片段在染色体上的位置， 每一列的可能又多，写到一半就懒了...

longbow0

可能会占用巨大内存的方法。具体消耗的内存与每条序列长度相关。

1. 
   
2. #!/usr/bin/perl
   
3. 
   
4. use strict;
   
5. use warnings;
   
6. 
   
7. use Bio::SeqIO;
   
8. 
   
9. my $seqi = Bio::SeqIO->new(
   
10.     -file => 'chr.fasta',
    
11.     -format => 'fasta',
    
12. );
    
13. 
    
14. while ( my $seq = $seqi->next_seq ) {
    
15.     my @sseqs = split /[a-z]+/, $seq->seq; # 提取大写序列
    
16.    
    
17.     # 或者更精确一点
    
18.     # my @sseqs = split /[atcg]+/, $seq->seq;
    
19. 
    
20.     for my $sseq ( @sseqs ) {
    
21.         next if length( $sseq ) < 500;
    
22. 
    
23.         print $sseq, "\n"; # 直接打印
    
24.     }
    
25. }
    
26. 
    

复制代码

mingxw

不耗内存版，行处理：

1. 
   
2. use strict;
   
3. use warnings;
   
4. 
   
5. my $repeat=$ARGV[0];#file repeat.fa
   
6. my $out=$ARGV[1];#splited file
   
7. my $min_seq_length=500;
   
8. my $out_line_length=50;
   
9. 
   
10. #repeat file load
    
11. open REPEAT,"<$repeat" or die "hi:$!";
    
12. open OUT,">$out" or die"Oh:$!";
    
13. my $DNA="";
    
14. my $ok_flag=1;
    
15. my $txt;
    
16. while($txt=<REPEAT>){
    
17.   chomp $txt;
    
18.   next unless($txt);
    
19.   if($txt =~ /^\>/){  #new DNA, clean the old one
    
20.     if($ok_flag) {
    
21.             while(length($DNA)>$out_line_length) {
    
22.                     print OUT substr($DNA,0,$out_line_length,""),"\n";
    
23.             }
    
24.             print OUT $DNA,"\n";
    
25.             $ok_flag=0;
    
26.     }
    
27.     elsif(length($DNA)>=$min_seq_length) {
    
28.             while(length($DNA)>$out_line_length) {
    
29.                     print OUT substr($DNA,0,$out_line_length,""),"\n";
    
30.             }
    
31.             print OUT $DNA,"\n";
    
32.     }
    
33.     $DNA="";
    
34.     print OUT "$txt\n";
    
35.   }
    
36.   elsif($txt=~/[atcg]/){
    
37.          
    
38.                 if($txt=~/^([ATCG]+)/) {
    
39.                         $DNA .=$1;
    
40.                 }
    
41.                
    
42.                 if($ok_flag) {
    
43.             while(length($DNA)>$out_line_length) {
    
44.                     print OUT substr($DNA,0,$out_line_length,""),"\n";
    
45.             }
    
46.             print OUT $DNA,"\n";       
    
47.             $ok_flag = 0;
    
48.     }
    
49.     else {
    
50.                         if(length($DNA)>=$min_seq_length) {
    
51.                     while(length($DNA)>$out_line_length) {
    
52.                             print OUT substr($DNA,0,$out_line_length,""),"\n";
    
53.                     }
    
54.                     print OUT $DNA,"\n";       
    
55.             }          
    
56.     }
    
57.                 $DNA="";
    
58.                 if($txt=~/([ATCG]+)$/) {
    
59.                         $DNA = $1;
    
60.           }         
    
61.         }
    
62.         else {
    
63.                 $DNA .= $txt;
    
64.                
    
65.                 if($ok_flag) {
    
66.             while(length($DNA)>$out_line_length) {
    
67.                     print OUT substr($DNA,0,$out_line_length,""),"\n";
    
68.             }
    
69.     }
    
70.     else {
    
71.                         if(length($DNA)>=$min_seq_length) {
    
72.                     while(length($DNA)>$out_line_length) {
    
73.                             print OUT substr($DNA,0,$out_line_length,""),"\n";
    
74.                     }
    
75.                     $ok_flag = 1;
    
76.             }          
    
77.     }               
    
78.         }
    
79. }
    
80. 
    
81. close REPEAT;
    
82. close OUT;
    
83. 
    

复制代码

chenhao392

回复 8# longbow0


    谢谢了， 跑了一遍，很快，1.6 G~2.4G内存之间，看来Seq:IO很强大啊..
    只是，我的scripts其实还计算了每一段大写片段的在基因组上的位置...用split ...不好记录吧：）