多行时两两之间的皮尔森相关系数

xingzhou823

下面是皮尔森相关系数的计算代码，是我copy的。它只是求两行的系数值的，但是我想求n行两两之间的相关系数。不知道怎么办了。

1. #!/usr/local/bin/perl -w
   
2. @x=(110.869429,224.432417,135.254959,314.799695,59.180641,160.191048,52.39372,283.666241,186.943458,170.444396,156.968469,424.092255,108.10017,51.204317,156.07176075,114.350089,66.6382395,154.04725,109.046664,228.484335
   
3. );
   
4. @y=(17.327247,17.404148,18.458939,17.559107,16.760933,15.49883,16.692633,17.327166,18.539049,19.794257,16.585954,15.774536,18.639974,16.643166,17.98409575,15.085853,18.899189,15.530345,16.620445,14.616374
   
5. );
   
6. if($#x==$#y){
   
7. $N = $#x;
   
8. print "N = $N\n";
   
9. $sum_sq_x = 0;
   
10. $sum_sq_y = 0;
    
11. $sum_coproduct = 0;
    
12. $mean_x = $x[1];
    
13. $mean_y = $y[1];
    
14. for($i=2;$i<=$N;$i++){
    
15.      $sweep = ($i - 1.0) / $i;
    
16.      $delta_x = $x[$i] - $mean_x;
    
17.      $delta_y = $y[$i] - $mean_y;
    
18.      $sum_sq_x += $delta_x * $delta_x * $sweep;
    
19.      $sum_sq_y += $delta_y * $delta_y * $sweep;
    
20.      $sum_coproduct += $delta_x * $delta_y * $sweep;
    
21.      $mean_x += $delta_x / $i;
    
22.      $mean_y += $delta_y / $i;
    
23. }
    
24. $pop_sd_x = sqrt( $sum_sq_x );
    
25. $pop_sd_y = sqrt( $sum_sq_y );
    
26. $cov_x_y = $sum_coproduct;
    
27. $correlation = $cov_x_y / ($pop_sd_x * $pop_sd_y);
    
28. print "Calculate Pearson correlation euqals:\t$correlation\n";
    
29. exit(0);
    
30. }
    
31. print "PASS!\n";
    
32. exit(1);
    

复制代码

wind_ch

n行兩兩之間的係數值 和 兩行之間的係數值有什麽關係

xingzhou823

回复 2# wind_ch


    我不知道呢。
我只想求n行中，所有每两行之间的系数值。
要是有10行，那么我就要求45个系数值。

wind_ch

和你程序里面一样, n行的话把每行保存成一个数组, 这是一个组合的问题吧.

xingzhou823

回复 4# wind_ch


可是，这个对我来说还是很难。
我卡在这儿了。
我就是得到两两之间的系数值，并输出
a   b   0.5
a   c    -0.3
b    d    0.8

wind_ch

假設你有10行數據, 每行保存進一個數組裏面, @arr1, @arr2 ... @arr10;
你可以每得到一個數組就把引用push到@data中去
把你處理兩行的程序寫成函數pearson

1. 
   
2. use strict;
   
3. 
   
4. sub pearson($){
   
5.     my ($ref_a, $ref_b) = @_;
   
6.     my @x = @{$ref_a};
   
7.     my @y = @{$ref_b};
   
8. 
   
9.     ...
   
10. }
    
11. 
    
12. my @data = (\@arr1, \@arr2, ... , \@arr10);
    
13. for(my $i = 0;  $i < $#data; $i++){
    
14.     for(my $j = $i + 1; $j <= $#data; $j++){
    
15.         pearson($data[$i], $data[$j]);
    
16.     }
    
17. }   
    

复制代码

xingzhou823

回复 6# wind_ch

1. #!/usr/local/bin/perl -w
   
2. 
   
3. sub pearson($){
   
4.         my ($ref_a,$ref_b) = @_;
   
5.         my @x = @{$ref_a};
   
6.         my @y = @{$ref_b};
   
7.         my $N;
   
8. if($#x==$#y){
   
9. $N = $#x;
   
10. print "N = $N\n";
    
11. my $sum_sq_x = 0;
    
12. my $sum_sq_y = 0;
    
13. my $sum_coproduct = 0;
    
14. my $mean_x = $x[1];
    
15. my $mean_y = $y[1];
    
16. for($i=2;$i<=$N;$i++){
    
17.      $sweep = ($i - 1.0) / $i;
    
18.      $delta_x = $x[$i] - $mean_x;
    
19.      $delta_y = $y[$i] - $mean_y;
    
20.      $sum_sq_x += $delta_x * $delta_x * $sweep;
    
21.      $sum_sq_y += $delta_y * $delta_y * $sweep;
    
22.      $sum_coproduct += $delta_x * $delta_y * $sweep;
    
23.      $mean_x += $delta_x / $i;
    
24.      $mean_y += $delta_y / $i;
    
25. }
    
26. $pop_sd_x = sqrt( $sum_sq_x );
    
27. $pop_sd_y = sqrt( $sum_sq_y );
    
28. $cov_x_y = $sum_coproduct;
    
29. $correlation = $cov_x_y / ($pop_sd_x * $pop_sd_y);
    
30. print "Calculate Pearson correlation euqals:\t$correlation\n";
    
31. exit(0);
    
32. }
    
33. print "PASS!\n";
    
34. exit(1);
    
35. }
    
36. 
    
37. open my $fh,"<","rece_test05.txt" or die "Error.\n";
    
38. my $out_file = $file . ".out";
    
39. open WFH,">",$out_file or die "Cannot open $out_file : $!\n";
    
40. 
    
41. 
    
42. my @data = (\@arr1, \@arr2, ... , \@arr500);
    
43. for(my $i = 0;  $i < $#data; $i++){
    
44.     for(my $j = $i + 1; $j <= $#data; $j++){
    
45.         pearson($data[$i], $data[$j]);
    
46.     }
    
47. }   

复制代码

我的文件是这样的，并且有500行，太多了。

1. LOC_Os04g46810.1         7.725284         7.728594         8.770031         9.910465         8.743908        10.332676
   
2. LOC_Os12g38010.1         50.067357        17.339407        23.806961        18.379239        22.397611        18.687081
   
3. LOC_Os10g31680.1         15.484767        14.557025        24.277390        22.325550        14.808402        21.870804
   
4. LOC_Os12g38310.1         5.636056         4.594795          5.666566         5.314397          5.830249        6.089081       
   
5. LOC_Os01g24560.1         7.853027         6.814309        11.343256         7.927583         9.649349        8.403233       
   
6. LOC_Os11g24140.1         11.791185        16.645000        16.399312        17.516789        29.046365        24.098895
   
7. LOC_Os10g31640.1         7.241884         6.794815         7.552480         8.329901         7.704593         11.929825
   
8. LOC_Os07g01410.1         20.902209        20.462937        26.375332        19.077440        23.761116        24.877955
   
9. LOC_Os10g18820.1          6.598861         8.858383         6.343684         7.583544         6.956356        8.908550       
   
10. LOC_Os01g24550.1         17.682395        18.022256        23.544243        18.772425        15.611792        20.254603
    

复制代码

我不知道，怎么把每一行引入数组@arr1，@arr2，，，而且，还得$name1 = pop(@arr1),$name2 = pop (arr2),得到每一行的 名称。
我想输出的文件格式是
LOC_Os12g38010.1    LOC_Os04g46810.1  0.356489
既 $name1 $name2 $correlation
关键是有太多行

wind_ch

這個你稍微看下書就能解決的

xingzhou823

回复 8# wind_ch


    这个也太难了吧，子程序的输出部分我不会，文件的读入我也不会，卡在这好难受啊

iamlimeng

1. #!/usr/bin/perl
   
2. 
   
3. use strict;
   
4. use warnings;
   
5. my $file = "rece_test05";
   
6. 
   
7. open (IN,"<$file.txt") or die "Error.\n";
   
8. my @data = <IN>;
   
9. close IN;
   
10. chomp @data;
    
11. my $out_file = $file . ".out";
    
12. open (OUT,">$out_file.txt") or die "Cannot open $out_file : $!\n";
    
13. my ($name1,$name2);
    
14. foreach my $i (0..$#data-1){
    
15.         foreach my $j($i+1..$#data) {
    
16.                  my @line1 = split(/\s+/,$data[$i]);
    
17.                  my @line2 = split(/\s+/,$data[$j]);
    
18.                  $name1 = shift @line1;
    
19.                  $name2 = shift @line2;
    
20.                  pearson(\@line1,\@line2);
    
21.          }
    
22. }
    
23. close OUT;
    
24. 
    
25. sub pearson {
    
26.     my ($ref_a, $ref_b) = @_;
    
27.     my @x = @{$ref_a};
    
28.     my @y = @{$ref_b};
    
29.         if($#x==$#y){
    
30.                  my $N = $#x;
    
31.                  my $sum_sq_x = 0;
    
32.                  my $sum_sq_y = 0;
    
33.                  my $sum_coproduct = 0;
    
34.                  my $mean_x = $x[1];
    
35.                  my $mean_y = $y[1];
    
36.                  for(my $i=2;$i<=$N;$i++){
    
37.                           my $sweep = ($i - 1.0) / $i;
    
38.                           my $delta_x = $x[$i] - $mean_x;
    
39.                           my $delta_y = $y[$i] - $mean_y;
    
40.                           $sum_sq_x += $delta_x * $delta_x * $sweep;
    
41.                           $sum_sq_y += $delta_y * $delta_y * $sweep;
    
42.                           $sum_coproduct += $delta_x * $delta_y * $sweep;
    
43.                           $mean_x += $delta_x / $i;
    
44.                           $mean_y += $delta_y / $i;
    
45.                  }
    
46.                  my $pop_sd_x = sqrt( $sum_sq_x );
    
47.                  my $pop_sd_y = sqrt( $sum_sq_y );
    
48.                  my $cov_x_y = $sum_coproduct;
    
49.                  my $correlation = $cov_x_y / ($pop_sd_x * $pop_sd_y);
    
50.                  print OUT "$name1 $name2: $correlation\n";
    
51.         }
    
52. }

复制代码