[Bio] / FigKernelPackages / ExpressionDir.pm Repository:
ViewVC logotype

Diff of /FigKernelPackages/ExpressionDir.pm

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.8, Mon Jan 24 19:42:26 2011 UTC revision 1.9, Fri Feb 4 22:55:11 2011 UTC
# Line 999  Line 999 
999      return $res;      return $res;
1000  }  }
1001    
1002  1;  sub best_pearson_corr {
1003        my($self,$pegs1,$cutoff) = @_;
1004    
1005        my @pegs2 = $self->all_features('peg');
1006        my $handle = $self->get_pc_hash_strip($pegs1,\@pegs2);
1007    
1008        my %ok;
1009        my $i;
1010        for ($i=0; ($i < @$pegs1); $i++)
1011        {
1012            foreach my $peg2 ( @pegs2 )
1013            {
1014                my $pc = &pearson_corr($handle,$pegs1->[$i],$peg2);
1015                if (abs($pc >= $cutoff))
1016                {
1017                    $ok{$pegs1->[$i]} -> {$peg2} = $pc;
1018                }
1019            }
1020        }
1021        return \%ok;
1022    }
1023    
1024    sub pearson_corr {
1025        my($hash,$peg1,$peg2) = @_;
1026        my $v = $hash->{$peg1}->{$peg2};
1027        return defined($v) ? sprintf("%0.3f",$v) : " ";
1028    }
1029    
1030    sub get_pc_hash_strip {
1031        my($self,$pegs1,$pegs2) = @_;
1032        my $corrH = $self->get_corr;
1033        my $hash  = &compute_pc_strip($pegs1,$pegs2,$corrH);
1034        return $hash;
1035    }
1036    
1037    sub get_corr {
1038        my($self) = @_;
1039    
1040        my $dir           = $self->expr_dir;
1041        my $rawF          = "$dir/rma_normalized.tab";
1042        my %gene_to_values;
1043        open(RAW,"<$rawF") || die "could not open $rawF";
1044        while (<RAW>)
1045        {
1046            chomp;
1047            my ($gene_id, @gxp_values) = split("\t");
1048            $gene_to_values{$gene_id} = \@gxp_values;
1049        }
1050        close(RAW);
1051        return \%gene_to_values;
1052    }
1053    
1054    sub compute_pc_strip {
1055        my ($pegs1,$pegs2, $gxp_hash) = @_;
1056        my %values = ();
1057    
1058        for (my $i = 0; $i < @$pegs1; $i++)
1059        {
1060            my $stat = Statistics::Descriptive::Full->new();
1061            $stat->add_data(@{$gxp_hash->{$pegs1->[$i]}});
1062    
1063            foreach my $peg2 (@$pegs2)
1064            {
1065                if ($pegs1->[$i] ne $peg2)
1066                {
1067                    my ($q, $m, $r, $err) = $stat->least_squares_fit(@{$gxp_hash->{$peg2}});
1068                    $values{$pegs1->[$i]}->{$peg2} = $r;
1069                }
1070            }
1071        }
1072    
1073        return \%values;
1074    }
1075    
1076    
1077    1;
1078    

Legend:
Removed from v.1.8  
changed lines
  Added in v.1.9

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3