对多个文件求平均

huang6894

还有一个问题。。。。对于N个以下的文本，我需要对第2、4、7列相同的第5列求平均值并作分析，想问一下应该怎么优化我的代码~

这N个文件的顺序格式都是一样的，不过每个文件在100M左右，每个文件里面的每一行都是唯一的。第2、第4列是位置信息，第3列则是区域信息，第四列是该区域的每个点，区域文件存放在文件ref上，第7列是正负信息，第5列则是对应的质量值。
设定第6列是质控值，用200除以该值得到权重，求平均时需要把第五列先乘以该权重才是他的真实质量值，每个文件的权重是唯一的。
目的是计算每个位置在所有文本中的平均值，然后分析假如某文件某行的数值小于0.5倍平均值的话，输出该行。。。。
文本：

1. 
   
2. ==> head -6 result_for_depth.13B0005773.1408345911.txt result_for_depth.13B0005805.1408342063.txt result_for_depth.13B0013031.1408342008.txt result_for_depth.13B0013194.1408342004.txt
   
3. 
   
4. ==> result_for_depth.13B0005773.1408345911.txt <==
   
5. 13B0005773      chr1    2336191-2337323 2336191 55      123     +
   
6. 13B0005773      chr1    2336191-2337323 2336191 16      123     -
   
7. 13B0005773      chr1    2336191-2337323 2336192 52      123     +
   
8. 13B0005773      chr1    2336191-2337323 2336192 16      123     -
   
9. 13B0005773      chr1    2336191-2337323 2336193 52      123     +
   
10. 13B0005773      chr1    2336191-2337323 2336193 16      123     -
    
11. 
    
12. ==> result_for_depth.13B0005805.1408342063.txt <==
    
13. 13B0005805      chr1    2336191-2337323 2336191 54      119     +
    
14. 13B0005805      chr1    2336191-2337323 2336191 11      119     -
    
15. 13B0005805      chr1    2336191-2337323 2336192 50      119     +
    
16. 13B0005805      chr1    2336191-2337323 2336192 11      119     -
    
17. 13B0005805      chr1    2336191-2337323 2336193 50      119     +
    
18. 13B0005805      chr1    2336191-2337323 2336193 11      119     -
    
19. 
    
20. ==> result_for_depth.13B0013031.1408342008.txt <==
    
21. 13B0013031      chr1    2336191-2337323 2336191 51      118     +
    
22. 13B0013031      chr1    2336191-2337323 2336191 14      118     -
    
23. 13B0013031      chr1    2336191-2337323 2336192 47      118     +
    
24. 13B0013031      chr1    2336191-2337323 2336192 14      118     -
    
25. 13B0013031      chr1    2336191-2337323 2336193 47      118     +
    
26. 13B0013031      chr1    2336191-2337323 2336193 15      118     -
    
27. 
    
28. ==> result_for_depth.13B0013194.1408342004.txt <==
    
29. 13B0013194      chr1    2336191-2337323 2336191 49      117     +
    
30. 13B0013194      chr1    2336191-2337323 2336191 13      117     -
    
31. 13B0013194      chr1    2336191-2337323 2336192 43      117     +
    
32. 13B0013194      chr1    2336191-2337323 2336192 13      117     -
    
33. 13B0013194      chr1    2336191-2337323 2336193 43      117     +
    
34. 13B0013194      chr1    2336191-2337323 2336193 13      117     -
    

复制代码

文件ref：

1. 
   
2. chr1    2336191    2337323
   

复制代码

我的代码：

1. 
   
2. die "-- must set depth_dir and ref_bed --\n" if (@ARGV < 2);
   
3. my $dir = $ARGV[0];
   
4. my $bed = $ARGV[1];
   
5. my(%dep,%sam,%f,%p);
   
6. my @allfile=glob "$ARGV[0]/result_for_depth.*.txt";
   
7. foreach my $file(@allfile){
   
8.         open FILE,$file;
   
9.         while(<FILE>){
   
10.                 chomp;
    
11.                  my ($sample, $chr, $pos, $dep, $qc, $stand) = (split /\t/,$_)[0,1,3,4,5,6];
    
12.                  my $depth = $dep * 200/$qc;
    
13.                  my $yu = $pos.$stand;
    
14.                  $dep{$chr}{$yu}+=$depth;
    
15.                  $sam{$file}{$chr}{$yu}=$depth;
    
16.                  $f{$file}{$chr}{$pos}="$sample\t$chr\t$pos";
    
17.          }
    
18.          close FILE;
    
19. }
    
20. 
    
21. foreach my $file(@allfile){
    
22.         open OUT,"> $file.cnv.result" || die "$!";
    
23.         open BED,"$bed" or die "$!";
    
24.         while(<BED>){
    
25.                 chomp;
    
26.                 my ($chr, $pos1, $pos2) = (split /\t/,$_)[0,1,2];
    
27.                 for my $pos ($pos1 .. $pos2){
    
28.                         my $y = $pos."+";
    
29.                         my $t = $pos."-";
    
30.                         my $mean1 = int(($dep{$chr}{$y}/$#allfile)+1);
    
31.                         my $mean2 = int(($dep{$chr}{$t}/$#allfile)+1);
    
32.                         if(($sam{$file}{$chr}{$y}<=0.5*$mean1)&&($sam{$file}{$chr}{$t}<=0.5*$mean2)){               
    
33.                                 print OUT "$f{$file}{$chr}{$pos}\n";
    
34.                         }
    
35.                   }
    
36.             }
    
37. }
    
38. 
    
39. 
    

复制代码

chenhao392

全读入hash么....要省内存，sort然后while(<>){}就好。

1. cat result_for_depth*txt | sort -k2 -k4n -k7 >working_file.bed

复制代码

然后

1. #!/usr/bin/perl
   
2. 
   
3. use strict;
   
4. use warnings;
   
5. 
   
6. 
   
7. my $pattern='';
   
8. my %hash;
   
9. my @f_score;
   
10. my @r_score;
    
11. 
    
12. open FILE, "<$ARGV[0]" or die "$. $!\n";
    
13. while(<FILE>){
    
14.     chomp;
    
15.     my ($sample, $chr, undef, $pos, $dep, $qc, $strand)=split(/\s+/,$_);
    
16.     if($pattern ne "$chr$pos"){
    
17.         $pattern = $chr.$pos;
    
18.         if($. >1){
    
19.             &flush();
    
20.         }
    
21.     }
    
22.     else{
    
23.     }
    
24. 
    
25.     my $score=$dep*200/$qc;
    
26.     if($strand eq '+'){
    
27.         push @f_score,$score;
    
28.         push @{$hash{'+'}},"$sample\t$chr\t$pos\t$score";
    
29.     }
    
30.     else{
    
31.         push @r_score,$score;
    
32.         push @{$hash{'-'}},"$sample\t$chr\t$pos\t$score";
    
33.     }
    
34. }
    
35. close FILE;
    
36. 
    
37. &flush();
    
38. 
    
39. 
    
40. 
    
41. 
    
42. sub flush{
    
43.     my $f_mean=mean(@f_score);
    
44.     my $r_mean=mean(@r_score);
    
45.     my %mem;
    
46.     foreach my $slide(@{$hash{'+'}}){
    
47.         my ($sample,$chr,$pos,$score)=split(/\t/,$slide);
    
48.         if($score >= 0.5*$f_mean){
    
49.             $mem{"$sample\t$chr\t$pos"}="";
    
50.         }
    
51.     }
    
52.     foreach my $slide(@{$hash{'-'}}){
    
53.         my ($sample,$chr,$pos,$score)=split(/\t/,$slide);
    
54.         if($score >= 0.5*$r_mean && defined $mem{"$sample\t$chr\t$pos"}){
    
55.             print "$sample\t$chr\t$pos\n";
    
56.         }
    
57.     }
    
58. 
    
59.     %mem=();
    
60.     @f_score=();
    
61.     @r_score=();
    
62.     %hash=();
    
63. }
    
64. 
    
65. sub mean{
    
66.     my $sum;
    
67.     foreach(@_){
    
68.         $sum+=$_;
    
69.     }
    
70.     return $sum/($#_+1);
    
71. }
    

复制代码

1. perl tt.pl working_file.bed
   
2. 13B0005805        chr1        2336191
   
3. 13B0013194        chr1        2336191
   
4. 13B0013031        chr1        2336191
   
5. 13B0005773        chr1        2336191
   
6. 13B0005805        chr1        2336192
   
7. 13B0013194        chr1        2336192
   
8. 13B0013031        chr1        2336192
   
9. 13B0005773        chr1        2336192
   
10. 13B0005805        chr1        2336193
    
11. 13B0013194        chr1        2336193
    
12. 13B0013031        chr1        2336193
    
13. 13B0005773        chr1        2336193
    

复制代码

pitonas

优化下 ~ {:2_172:}
            my $mean2 = int( ( $dep{$chr}{$t} / $#allfile ) + 1 );

$#allfile: index ==> count - 1
@allfile: count

1. #!/usr/bin/perl
   
2. 
   
3. die "-- must set depth_dir and ref_bed --\n" if ( @ARGV < 2 );
   
4. my $dir = $ARGV[0];
   
5. my $bed = $ARGV[1];
   
6. my ( %dep, %sam, %f, %p );
   
7. my @allfile = glob "$ARGV[0]/result_for_depth.*.txt";
   
8. 
   
9. foreach my $file (@allfile) {
   
10.     open FILE, $file;
    
11.     while (<FILE>) {
    
12.         my ( $sample, $chr, undef, $pos, $dep, $qc, $stand ) = split;
    
13.         my $depth = $dep * 200 / $qc;
    
14.         my $yu    = $pos . $stand;
    
15.         $dep{$chr}{$yu} += $depth;
    
16.         $sam{$file}{$chr}{$yu} = $depth;
    
17.         $f{$file}{$chr}{$pos}  = "$sample\t$chr\t$pos";
    
18.     }
    
19.     close FILE;
    
20. }
    
21. 
    
22. open BED, "$bed" or die "$!";
    
23. my @BED = map [split], <BED>;
    
24. 
    
25. foreach my $file (@allfile) {
    
26.     open OUT, "> $file.cnv.result" || die "$!";
    
27.     for (@BED) {
    
28.         my ( $chr, $pos1, $pos2 ) = @$_;
    
29.         for my $pos ( $pos1 .. $pos2 ) {
    
30.             my $y     = $pos . "+";
    
31.             my $t     = $pos . "-";
    
32.             
    
33.             # $#allfile: index, @allfile: count
    
34.             my $mean1 = int( ( $dep{$chr}{$y} / $#allfile ) + 1 );
    
35.             my $mean2 = int( ( $dep{$chr}{$t} / $#allfile ) + 1 );
    
36.             if (   ( $sam{$file}{$chr}{$y} <= 0.5 * $mean1 )
    
37.                 && ( $sam{$file}{$chr}{$t} <= 0.5 * $mean2 ) )
    
38.             {
    
39.                 print OUT "$f{$file}{$chr}{$pos}\n";
    
40.             }
    
41.         }
    
42.     }
    
43. }
    
44. 
    

复制代码

huang6894

回复 2# pitonas


    嗯嗯，谢谢pitonas大大，我研究一下

huang6894

回复 2# pitonas


    大神。。。程序死掉了，6G内存跑不了。。。。

pitonas

我不知道它是否可以改写?
我是文盲, 看不懂问题. 我做不到了....


非常抱歉 ~ {:2_168:}

回复 4# huang6894


   

huang6894

回复 5# pitonas


    好吧，我又错了。。。。。

我修改了一下：
不过还在测试。。。希望有好运气吧。。。

1. die "-- must set depth_dir and ref_bed --\n" unless (-d $ARGV[0] );
   
2. my $dir = $ARGV[0];
   
3. my %dep;
   
4. my @allfile = glob "$dir/result_for_depth.*.txt";
   
5. 
   
6. foreach my $file (@allfile) {
   
7.         my $line;
   
8.         open FILE, $file;
   
9.         while (<FILE>) {
   
10.                 chomp;
    
11.                 my ( $sample, $chr, undef, $pos, $depth, $stand ) = split;#输入文件修改了，在原有的基础上做了权重校正，这里去掉了$qc
    
12.                 push @{$dep{$line}{R}},[$sample,$depth];#因为所有文件的每一行的key都是一样的
    
13.                 $dep{$line}{$sample}=[ $sample, $chr, $pos, $depth, $stand ];
    
14.                 $line ++;
    
15.                }
    
16.                close FILE;
    
17. }
    
18. 
    
19. open OUT,"> $dir/all_cnv_result";
    
20. while ( my ( $k, $v ) = each %dep ) {
    
21.         my @tmp = @{ $v->{R} };
    
22.         my ($a,@R);
    
23.         for my $t (@tmp) {
    
24.                  $a += $t->[1];
    
25.                  push @R,$t->[0];
    
26.         }
    
27.         my $mean = $a/($#allfile+1);
    
28.         for my $R (@R) {
    
29.                 my ( $A, $B, $C, $D, $E ) = @{ $v->{ $R } };
    
30.                 if($D <= 0.45*$mean){
    
31.                         print OUT "$A\t$B\t$C\t$D\t$E\t$mean\tdel\n";
    
32.                 }
    
33.         }
    
34. }
    
35. close OUT;
    

复制代码

pitonas

1.                 # chomp; # don't do this !! it's stupid!!
   
2.                 my ( $sample, $chr, undef, $pos, $depth, $stand ) = split;

复制代码

回复 6# huang6894


   

huang6894

回复 7# pitonas


    哦。。。好的

yestreenstars

回复 7# pitonas

why stupid?