perl的內(nèi)存占用，求教&探討

chenhao392

今天有個生物信息的任務(wù)，完成了，但是很不爽的搞死了電腦兩次。于是全程播報我的思路和步驟，來論壇求教了，求拍磚：

1. 數(shù)據(jù)描述

有基因組數(shù)據(jù)如下，大寫的是正常序列，小寫的是repeat區(qū)域：

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. 任務(wù)要求
按照小寫的序列區(qū)域，將一個長的染色體分割成較短的片段，小寫區(qū)域丟棄，大寫片段短于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. 測試小量數(shù)據(jù)通過
......


6. 大數(shù)據(jù)量，死機
于是在應(yīng)用在了基因組上，數(shù)據(jù)量2G多

瞬間占用40G以上內(nèi)存，服務(wù)器死機

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的內(nèi)存，程序順利跑完，但是我依然不明白哪里需要4G內(nèi)存了。
在我看來，這兩行：
$seq="";
undef(@line);
已經(jīng)釋放了內(nèi)存，應(yīng)該只是%chrom存儲的2G序列而已了。

zhlong8

2G 的原始數(shù)據(jù)才占 4G 的內(nèi)存，已經(jīng)很省了啊？

chenhao392

回復 2# zhlong8


    這個....這樣啊....
    那另外2G內(nèi)存被什么吃掉了？

    或者說，
    為什么4G就很省了？

zhlong8

將一個長的染色體分割成較短的片段，小寫區(qū)域丟棄，大寫片段短于500字符就丟棄。

這個用 RE 就能做了��？你為什么一個一個字符處理？

@seq 那里 split 成單個字符每1K個占用內(nèi)存是 48K，如果你一段基因超長的話……

Perl 的字符串可變，所以可能分配的內(nèi)存比需要的多，而且是以2^n 的速度擴大的，所以平均長度越長浪費的也就越多

數(shù)據(jù)結(jié)構(gòu)不也要占用內(nèi)存，建一個空空的數(shù)組或 hash 難道就不需要空間了？

chenhao392

回復 4# zhlong8


謝謝說明！


將一個長的染色體分割成較短的片段，小寫區(qū)域丟棄，大寫片段短于500字符就丟棄。
這個用 RE 就能做了�。磕銥槭裁匆粋�一個字符處理？

我只能說是，想到什么就寫什么了，另，RE是什么？


@seq 那里 split 成單個字符每1K個占用內(nèi)存是 48K，如果你一段基因超長的話……

是的，恐怖的長，2G只有10個基因組，平均一個200MB. 看來這就是我第一次程序掛掉的原因


Perl 的字符串可變，所以可能分配的內(nèi)存比需要的多，而且是以2^n 的速度擴大的，所以平均長度越長浪費的也就越多
數(shù)據(jù)結(jié)構(gòu)不也要占用內(nèi)存，建一個空空的數(shù)組或 hash 難道就不需要空間了？

$seq的長度經(jīng)常會有幾千上萬，一直在變，那這就是原因了？
我是不是可以理解為：

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


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

前者比后者 消耗的內(nèi)存多？

zhlong8

對于很長的數(shù)據(jù)按行存儲和處理比較高效。

說用 RE 指的是對于每一行的數(shù)據(jù)可以用正則表達式來判斷分割。

對于比較短的字符串不需要考慮這么多，一個變量保存字符串的長度是曾經(jīng)最長的長度，即使賦值了個比較短的字符串，相當于剩余的空間就是緩沖區(qū)了。用 undef 可以釋放

chenhao392

回復 6# zhlong8


謝謝，變量占用內(nèi)存大概明白了...


    RE指，類似這樣的處理么？
     if (   $line=~/([ATCG]*)([atcg]*)([ATCG]*)/ ) {
                $1;
                $2;
                $3;
   }

其實一開始我是有想這么寫的，但是還需要計算大寫片段在染色體上的位置， 每一列的可能又多，寫到一半就懶了...

longbow0

可能會占用巨大內(nèi)存的方法。具體消耗的內(nèi)存與每條序列長度相關(guān)。

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

不耗內(nèi)存版，行處理：

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內(nèi)存之間，看來Seq:IO很強大啊..
    只是，我的scripts其實還計算了每一段大寫片段的在基因組上的位置...用split ...不好記錄吧：）