[Bio] / FigKernelScripts / ma.cgi Repository:
ViewVC logotype

View of /FigKernelScripts/ma.cgi

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (download) (annotate)
Tue Feb 9 02:43:31 2010 UTC (10 years ago) by overbeek
Branch: MAIN
CVS Tags: mgrast_dev_08112011, mgrast_dev_08022011, myrast_rel40, mgrast_dev_05262011, mgrast_dev_04082011, rast_rel_2010_0928, mgrast_version_3_2, mgrast_dev_12152011, mgrast_dev_06072011, rast_rel_2010_0526, mgrast_dev_10262011, mgrast_dev_02212011, rast_rel_2010_1206, mgrast_release_3_0, mgrast_dev_03252011, rast_rel_2011_0119, mgrast_release_3_0_4, mgrast_release_3_0_2, mgrast_release_3_0_3, mgrast_release_3_0_1, mgrast_dev_03312011, mgrast_release_3_1_2, mgrast_release_3_1_1, mgrast_release_3_1_0, mgrast_dev_04132011, mgrast_dev_04012011, rast_rel_2010_0827, myrast_33, rast_rel_2011_0928, mgrast_dev_04052011, mgrast_dev_02222011
added massaging and voting

#########################################################################
# -*- perl -*-
#
# Copyright (c) 2003-2006 University of Chicago and Fellowship
# for Interpretations of Genomes. All Rights Reserved.
#
# This file is part of the SEED Toolkit.
# 
# The SEED Toolkit is free software. You can redistribute
# it and/or modify it under the terms of the SEED Toolkit
# Public License. 
#
# You should have received a copy of the SEED Toolkit Public License
# along with this program; if not write to the University of Chicago
# at info@ci.uchicago.edu or the Fellowship for Interpretation of
# Genomes at veronika@thefig.info or download a copy from
# http://www.theseed.org/LICENSE.TXT.
#

use FIG;
my $fig = new FIG;

use DB_File;
use SSserver;
use SeedEnv;

use HTML;
use strict;

use CGI;
my $cgi = new CGI;

if (0)
{
    my $VAR1;
    eval(join("",`cat $FIG_Config::temp/ma_cgi`));
    $cgi = $VAR1;
    print STDERR &Dumper($cgi);
}

if (0)
{
    print $cgi->header;
    my @params = $cgi->param;
    print "<pre>\n";
    foreach $_ (@params)
    {
	print "$_\t:",join(",",$cgi->param($_)),":\n";
    }

    if (0)
    {
	if (open(TMP,">$FIG_Config::temp/ma_cgi"))
	{
	    print TMP &Dumper($cgi);
	    close(TMP);
	}
    }
    exit;
}
my $html = [];
my $source_dir = "/home/fig/Ross/MA/CGI";
my $dir        = "$FIG_Config::data/Regulons/Shewanella";
my $pearson_split = $cgi->param('pearson_split');
if (! defined($pearson_split)) { $pearson_split = 0.4 }
if (! -d $dir)
{
    &build_regulon_dir($source_dir,$dir,$pearson_split);
}

my $request       = $cgi->param('request');
my $cond          = $cgi->param('cond');
my $bad           = $cgi->param('bad');
my $pegs          = $cgi->param('pegs');
my $resetG        = $cgi->param('resetG');

if ($request && ($request eq 'reset'))
{
    &build_regulon_dir($source_dir,$dir,$pearson_split);
}

if ($request eq "show_conditions")
{
    push(@$html,`cat $dir/conditions.html`);
}
elsif ($pegs)
{
    &show_pegs($cgi,$html,$pegs,$dir);
}
elsif (($request eq "show_condition") && defined($cond))
{
    my($condition,$nameC,$desc) = map { (($_ =~ /^(\d+)\s+(\S+)\s+(\S.*\S)/) && (($cond eq $1) || $cond eq $2)) ? 
					    ($1,$2,$3) : () 
				      } `cat $dir/condition.index`;
    if ($nameC)
    {
	push(@$html,$cgi->h3("$condition: $nameC - $desc"),"<br><br>",&link_to_conditions);
    }
}
elsif (($request eq "show_bad") && $bad)
{
    &show_bad_regulon($fig,$cgi,$html,$bad,$dir);
}
elsif (($request eq "split_one") && $bad)
{
    &split_one_regulon($fig,$cgi,$html,$bad,$dir);
}
else
{
    my($vecs,undef) = &read_vecs($dir);
    my $bad = &get_bad_regulons($dir,$vecs,$cgi);
    &display_bad($fig,$cgi,$html,$bad,$dir);
}
&HTML::show_page($cgi,$html);
exit;

sub link_to_conditions {
    return "<a href=./ma.cgi?request=show_conditions>show condition table</a>";
}

sub build_regulon_dir {
    my($source_dir,$dir,$pearson_split) = @_;

    if (-d $dir) { system "rm -r $dir" }
    &FIG::verify_dir($dir);
    system "cp $source_dir/predicted.operons $source_dir/*.assertions $source_dir/condition.index $source_dir/conditions.html $source_dir/ma $source_dir/peg.corr $dir";
    system "pushd $dir > /dev/null; ln -s /home/fig/Ross/MA/CGI/pearson.DB; popd > /dev/null";
    system "chmod -R 777 $dir";

    &FIG::verify_dir("$dir/Regulons");
    opendir(DIR,$dir) || die "could not open $dir";
    my @assF = grep { $_ =~ /\.assertions$/ } readdir(DIR);
    closedir(DIR);
    my $n = 1;
    foreach my $file (@assF)
    {
	open(IN,"<$dir/$file") || die "could not open $file";
	while (defined($_ = <IN>))
	{
	    if ($_ =~ /^(\S+)\t(\S.*\S)/)
	    {
		my $pegs = $1;
		my $rest = $2;
		my %hash = map { $_ => 1 } split(/,/,$pegs);
		my @in = sort { &FIG::by_fig_id($a,$b) } keys(%hash);
		if (@in > 1)
		{
		    open(OUT,">$dir/Regulons/$n") || die "could not open $dir/Regulons/$n";
		    my @split = &split_by_pearson($dir,\@in,$pearson_split);
		    foreach my $set (@split)
		    {
			print OUT join(",",@$set),"\t$rest\n";
		    }
		    close(OUT);
		    $n++;
		}
	    }
	}
	close(IN);
    }
}

sub split_by_pearson {
    my($dir,$pegs,$pearson_split) = @_;

#    print STDERR &Dumper($pegs);
    my($pearsonH,$pearson_hash_tieP) = &get_pearson_DB($dir);
    my %conn;
    my($i,$j,$peg1,$peg2,$p1,$p2);
    for ($i=0; ($i < @$pegs); $i++)
    {
	$peg1 = $pegs->[$i];
	$peg1 =~ /(\d+)$/;
	$p1 = $1;

	for ($j=$i+1; ($j < @$pegs); $j++)
	{
	    if ($i != $j)
	    {
		$peg2 = $pegs->[$j];
		$peg2 =~ /(\d+)$/;
		$p2 = $1;
		my $key = ($p1 < $p2) ? "$p1,$p2" : "$p2,$p1";
		my $x = $pearsonH->{$key};
		if (defined($x) && ($x >= $pearson_split))
		{
		    push(@{$conn{$peg1}},$peg2);
		    push(@{$conn{$peg2}},$peg1);
		}
	    }
	}
    }
#    print STDERR &Dumper(\%conn);

    my %used;
    my @sets = ();
    my @to_check = @$pegs;
    while (my $peg = shift @to_check)
    {
	next if ($used{$peg});

	my %in;
	my $set = [$peg];
	$used{$peg} = 1;
	$in{$peg}   = 1;
	for ($i=0; ($i < @$set); $i++)
	{
	    my $x = $conn{$set->[$i]};
	    foreach my $peg2 (@$x)
	    {
		if (! $in{$peg2})
		{
		    push(@$set,$peg2);
		    $in{$peg2} = 1;
		    $used{$peg2} = 1;
		}
	    }
	}
	push(@sets,$set);
    }
    
    &untie_pearson_hash($pearsonH,$pearson_hash_tieP);
#    print STDERR &Dumper(\@sets); die "aborted";
    return @sets;
}

sub get_pearson_DB {
    my($dir) = @_;
    my $pearson_hash;
    my $pearson_hash_tie = tie %$pearson_hash,'DB_File',"$dir/pearson.DB",O_RDONLY,0666,$DB_BTREE;
    $pearson_hash_tie || die "tie failed";
    return ($pearson_hash,\$pearson_hash_tie);
}

sub untie_pearson_hash {
    my($pearson_hash,$pearson_hash_tieP) = @_;

    undef $$pearson_hash_tieP;
    untie %$pearson_hash;
}

sub read_vecs {
    my($dir) = @_;

    my %seen;
    my $vecs = {};
    my $sz;
    my %val = ( '0' => 1, '1' => 2, '-1' => 0 );
    foreach my $x (map { $_ =~ /^([^\t]*\t){2}(fig\|\d+\.\d+\.peg\.\d+)([^\t]*\t){3}(\S.*\S)/; 
			 [$2,[map { $val{$_} } split(/\t/,$4)]] } 
		   `cat $dir/ma`)
    {
	if (! $sz)
	{
	    $sz = @{$x->[1]};
	}

	if (! $seen{$x->[0]})
	{
	    $seen{$x->[0]} = 1;
	    if (@{$x->[1]} != $sz)  { die 'malformed vectors' }
	    $vecs->{$x->[0]} = $x->[1];
	}
    }
    return ($vecs,$sz);
}

sub get_bad_regulons {
    my($dir,$vecs,$cgi) = @_;

    my $minsc = $cgi->param('minsc');
    if (! defined($minsc)) { $minsc = 6 }
    my $minON = $cgi->param('minON');
    if (! defined($minON)) { $minON = 0 }
    my $minOFF = $cgi->param('minOFF');
    if (! defined($minOFF)) { $minOFF = 0 }

    my $peg = $cgi->param('peg');
    my $bad = {};
    opendir(REGULONS,"$dir/Regulons") || die "could not open $dir/REGULONS";
    my @entries = grep { $_ =~ /^\d+$/ } readdir(REGULONS);
    closedir(REGULONS);
    foreach my $file (@entries)
    {
	open(IN,"<$dir/Regulons/$file") || die "could not open $dir/Regulons/$file";
	my @pegs = ();
	while (defined($_ = <IN>))
	{
	    if ($_ =~ /^(\S+)/)
	    {
		push(@pegs, split(/,/,$1));
	    }
	}
	close(IN);

	my($i,$j);
	my @pegs = sort { &FIG::by_fig_id($a,$b) } @pegs;

	for ($i=0; ($i < @pegs) && 
	           ((! $minON)  || (! $vecs->{$pegs[$i]}) || (@{&expressed_in($vecs->{$pegs[$i]})} >= $minON)) &&
	           ((! $minOFF) || (! $vecs->{$pegs[$i]}) || (@{&not_expressed_in($vecs->{$pegs[$i]})} >= $minOFF));
	     $i++) {}
	if ($i == @pegs)
	{
	    if ($peg && (! &member($peg,\@pegs))) { next }
	
	    for ($i=0; ($i < $#pegs); $i++)
	    {
		for ($j=$i+1; ($j < @pegs); $j++)
		{
		    my($v1,$v2);
		    if (($v1 = $vecs->{$pegs[$i]}) &&
			($v2 = $vecs->{$pegs[$j]}))
		    {
			my $sc = &scored($v1,$v2);
			if ($sc >= $minsc)
			{
			    $bad->{$file} = $sc;
			}
		    }
		}
	    }
	}
    }
    return $bad;
}

sub member {
    my($x,$xL) = @_;

    my $i;
    for ($i=0; ($i < @$xL) && ($x ne $xL->[$i]); $i++) {}
    return ($i < @$xL);
}

sub scored {
    my($v1,$v2) = @_;

    my $i;
    my $match = 0;
    my $mismatch = 0;
    my $sc = 0;
    for ($i=0; ($i < @$v1); $i++)
    {
	if (($v1->[$i] != 1) && ($v2->[$i] != 1))
	{
	    if ($v1->[$i] != $v2->[$i])
	    {
		$mismatch++;
	    }
	    else
	    {
		$match++;
	    }
	}
    }
    return ($match+$mismatch) ? int(($mismatch * 100) / ($match+$mismatch)) : 0;
}
	
sub display_bad {
    my($fig,$cgi,$html,$bad,$dir) = @_;

    my $n = 1;
    my $col_hdrs = ['','sc','regulon-id','Source'];
    my $tab = [];
    foreach my $file (sort { $bad->{$b} <=> $bad->{$a} } keys(%$bad))
    {
	open(TMP,"<$dir/Regulons/$file") || die "could not open $dir/Regulons/$file";
	$_ = <TMP>;
	chop;
	my @flds = split(/\t/,$_);
	shift @flds;
	my $extra = join("; ",@flds);
	close(TMP);
	push(@$tab,[$n,$bad->{$file},"<a href=./ma.cgi?request=show_bad&bad=$file>$file</a>",$extra]);
	$n++;
    }
    push(@$html,&HTML::make_table($col_hdrs,$tab,"Regulons that Need to be Checked"));
}

sub show_bad_regulon {
    my($fig,$cgi,$html,$bad,$dir) = @_;

    my $red = $cgi->param('red');
    if (! defined($red)) { $red = 0.001 }

    my $green = $cgi->param('green');
    if (! defined($green)) { $green = 1 }

    my($vecs,$sz) = &read_vecs($dir);

    my @links_to_cond;
    my $i;
    for ($i=0; ($i < $sz); $i++)
    {
	push(@links_to_cond,"<a target=condition href=ma.cgi?request=show_condition&cond=$i>*</a>");
    }

    my @peg_sets = map { [$_ =~ /fig\|\d+\.\d+\.peg\.\d+/g] } `cat $dir/Regulons/$bad`;
    my @all      = map { @$_ } @peg_sets;
    &show_correlation($html,\@all,$vecs);
    push(@$html,"<hr>");
    &show_pearson_correlation($dir,$html,\@all);
    &show_bad_orig($fig,$cgi,$html,$bad,$dir,$red,$green,$vecs,$sz,\@links_to_cond,\@all,"Putative Regulon");
    push(@$html,$cgi->hr,$cgi->h1('Data After Operon Analysis'));

    foreach my $peg_set (@peg_sets)
    {
	my $vecs1 = &massage_exp_vecs($fig,$vecs,$peg_set,$dir);
	&show_bad_orig($fig,$cgi,$html,$bad,$dir,$red,$green,$vecs1,$sz,\@links_to_cond,$peg_set,"Massaged Regulon");
    }
    push(@$html,$cgi->hr,$cgi->h1('Data After Voting'));

    foreach my $peg_set (@peg_sets)
    {
        my $vecs1 = &vote($vecs,$peg_set);
	&show_bad_orig($fig,$cgi,$html,$bad,$dir,$red,$green,$vecs1,$sz,\@links_to_cond,$peg_set,"After Voting");
    }
}

sub vote {
    my($vecs,$pegs) = @_;

    my %vecs1 = map { $vecs->{$_} ? ($_ => $vecs->{$_}) : () } @$pegs;
    my @ok_pegs = keys(%vecs1);
    if (@ok_pegs < 2) { return \%vecs1 }
    my $i;
    my $sz = @{$vecs1{$ok_pegs[0]}};
   for ($i=0; ($i < $sz); $i++)
    {
	my $c0 = 0;
	my $c2 = 0;
	foreach my $peg (@ok_pegs)
	{
	    my $c = $vecs1{$peg}->[$i];
	    if    ($c == 0) { $c0++ }
	    elsif ($c == 2) { $c2++ }
	}

	if ($c0 && ($c0 > $c2))
	{
	    &set_col(\%vecs1,$i,0);
	}

	if ($c2 && ($c0 < $c2))
	{
	    &set_col(\%vecs1,$i,2);
	}
    }
    return \%vecs1;
}

sub set_col {
    my($vecs1,$i,$v) = @_;

    foreach my $peg (keys(%$vecs1))
    {
	$vecs1->{$peg}->[$i] = $v;
    }
}

sub massage_exp_vecs {
    my($fig,$vecs,$pegs,$dir) = @_;
    
    my %vecs1 = map { $vecs->{$_} ? ($_ => $vecs->{$_}) : () } @$pegs;
    my @ok_pegs = keys(%vecs1);
    my @runs  = &get_runs($fig,\@ok_pegs,$dir,\%vecs1);
    &mark_changes_due_to_runs($fig,\@runs,\%vecs1);
    
    return \%vecs1;
}

sub get_runs {
    my($fig,$pegs,$dir,$vecs1) = @_;

    my %ops;
    my %relevant = map { $_ => 1 } @$pegs;
    my %sets     = map { $_ =~ /^(\d+)\t(\S+)/; ($1 && $relevant{$2}) ? ($1 => 1) : () }
                   `cat $dir/predicted.operons`;
    foreach $_ (`cat  $dir/predicted.operons`)
    {
	if (($_ =~ /^(\d+)\t(\S+)/) && $sets{$1} && $vecs1->{$2})
	{
	    push(@{$ops{$1}},$2);
	}
    }
    my @runs = map { my $x = $ops{$_}; (@$x > 1) ? [&sort_by_dir($fig,$x)] : () } keys(%ops);
    return @runs;
}

sub sort_by_dir {
    my($fig,$pegs) = @_;

    my @with_loc = map { my(undef,$beg,$end) = $fig->boundaries_of(scalar $fig->feature_location($_));
			 [$beg,$end,$_] } @$pegs;
    my @sorted   = sort { ($a->[0] < $a->[1]) ? ($a->[0] <=> $b->[0]) : ($b->[0] <=> $a->[0]) } @with_loc;
    return map { $_->[2] } @sorted;
}

sub mark_changes_due_to_runs {
    my($fig,$runs,$vecs1) = @_;

    foreach my $run (@$runs)
    {
	my $sz = @{$vecs1->{$run->[0]}};
	my $i;
	for ($i=0; ($i < $sz); $i++)
	{
	    &process_col($run,$vecs1,$i);
	}
    }
    return;
}

sub process_col {
    my($run,$vecs1,$i) = @_;

    my $j;
    for ($j=0; ($j < @$run) && ($vecs1->{$run->[$j]}->[$i] == 1); $j++) {}
    if ($j < @$run)
    {
	my $k;
	for ($k=0; ($k < @$run); $k++)
	{
	    $vecs1->{$run->[$k]}->[$i] = $vecs1->{$run->[$j]}->[$i];
	}
    }
}

sub show_bad_orig {
    my($fig,$cgi,$html,$bad,$dir,$red,$green,$vecs,$sz,$links_to_cond,$pegs,$hdr) = @_;
    my $col_hdrs = ['In New Set','PEG','aliases','Function',@$links_to_cond];
    my $tab = [];
    push(@$html, $cgi->start_form(-action => "ma.cgi"),
	 $cgi->hidden(-name => 'request',  -value => 'split_one', -override => 1),
	 $cgi->hidden(-name => 'bad',  -value => $bad));
    
    my @pvecs = map { $vecs->{$_} ? $vecs->{$_} : () } @$pegs;
    my $color = {};
    my($i,$j);
    for ($i=0; ($i < $sz); $i++)
    {
	my $on = 0;
	my $off = 0;
	for ($j=0; ($j < @pvecs); $j++)
	{
	    my $c = $pvecs[$j]->[$i];
	    if ($c == 0) { $off++ }
	    if ($c == 2) { $on++ }
	}
	if ($off && $on && 
	    (&FIG::min($off,$on) / ($off + $on) >= $red))
	{
	    $color->{$i} = 'red';
	}
	elsif (($off || $on) && (($off+$on) > 1) &&
	       ((&FIG::max($off,$on) / ($off + $on)) >= $green))
	{
	    $color->{$i} = 'green';
	}
    }
    foreach my $peg (@$pegs)
    {
	my $func = $fig->function_of($peg);
	my $v = $vecs->{$peg};
	if (! $v)
	{
	    $v = [];
	    for ($i=0; ($i < $sz); $i++)
	    {
		push(@$v,"");
	    }
	}
	my $checkbox = $cgi->checkbox( -name => 'checked', -value => $peg, -override => 1);
	my $aliases = $fig->feature_aliases($peg);
	my $base = [$checkbox,&peg_link($peg),$aliases,$func];
	
	for ($i=0; ($i < $sz); $i++)
	{
	    if ($color->{$i})
	    {
		push(@$base,[$v->[$i],"td bgcolor=$color->{$i}"]);
	    }
	    else
	    {
		push(@$base,$v->[$i]);
	    }
	}
	push(@$tab,$base);
    }
    push(@$html,&HTML::make_table($col_hdrs,$tab,$hdr));
    push(@$html,$cgi->submit('Split This One'),
	        $cgi->end_form);
}

sub peg_link { 
    my($peg) = @_;

    my $window = "PEG.$$";
    return "<a target=$window href=seedviewer.cgi?page=Annotation&feature=$peg>$peg</a>";
}

sub split_one_regulon {
    my($fig,$cgi,$html,$bad,$dir) = @_;

    my @checked = $cgi->param('checked');
    if (@checked == 0)
    {
	push(@$html,$cgi->h3("You need to check some pegs to achieve a split"));
    }
    else
    {
	my @old = `cat $dir/Regulons/$bad`;
	my $new;
	if ($new = &fix_reg(\@checked,\@old))
	{
	    if (@$new > 0)
	    {
		open(NEW,">$dir/Regulons/$bad") || die "could not overwrite $dir/Regulons/$bad";
		foreach $_ (@$new)
		{
		    print NEW $_;
		}
		close(NEW);
		push(@$html,$cgi->h3("Split $bad"));
	    }
	    else
	    {
		unlink("$dir/Regulons/$bad");
		push(@$html,$cgi->h3("deleted putative regulon set $bad"));
	    }
	}
	else
	{
	    push(@$html,$cgi->h3("split of $bad failed"));
	}
    }
}

sub fix_reg {
    my($checked,$old) = @_;

    my $new = [];
    my $seen_checked = 0;
    my %checked = map { $_ => 1 } @$checked;
    foreach my $line (@$old)
    {
	chop $line;
	my @flds = split(/\t/,$line);
	my $pegs = shift @flds;
	my @pegs = split(/,/,$pegs);
	my @s1 = grep { $checked{$_} } @pegs;
	my @s2 = grep { ! $checked{$_} } @pegs;
	if (@s2 > 0) { $seen_checked = 1 }
	if (@s2 > 1)
	{
	    push(@$new, join("\t",(join(",",@s2),@flds)) . "\n");
	}

	if (@s1 > 1)
	{
	    push(@$new, join("\t",(join(",",@s1),@flds)) . "\n");
	}
    }
    return $seen_checked ? $new : undef;
}
    
sub show_pegs {
    my($cgi,$html,$pegs,$dir) = @_;

    my @pegs = split(/,/,$pegs);
    my($vecs,$sz) = &read_vecs($dir);
    &show_correlation($html,\@pegs,$vecs);
    my($i,$j);
    my @links_to_cond;
    for ($i=0; ($i < $sz); $i++)
    {
	push(@links_to_cond,"<a target=condition href=ma.cgi?request=show_condition&cond=$i>*</a>");
    }
    my $col_hdrs = ['PEG','Function',@links_to_cond];
    my $tab = [];
    my @pvecs = map { $vecs->{$_} ? $vecs->{$_} : () } @pegs;
    my $color = {};
    for ($i=0; ($i < $sz); $i++)
    {
	my $on = 0;
	my $off = 0;
	for ($j=0; ($j < @pvecs); $j++)
	{
	    my $c = $pvecs[$j]->[$i];
	    if ($c == 0) { $off++ }
	    if ($c == 2) { $on++ }
	}
	if ($off && $on)
	{
	    $color->{$i} = 1;
	}
    }
    foreach my $peg (@pegs)
    {
	my $func = $fig->function_of($peg);
	my $v = $vecs->{$peg};
	if (! $v)
	{
	    $v = [];
	    for ($i=0; ($i < $sz); $i++)
	    {
		push(@$v,"");
	    }
	}
	my $base = [&peg_link($peg),$func];
	for ($i=0; ($i < $sz); $i++)
	{
	    if ($color->{$i})
	    {
		push(@$base,[$v->[$i],"td bgcolor=red"]);
	    }
	    else
	    {
		push(@$base,$v->[$i]);
	    }
	}
	push(@$tab,$base);
    }
    push(@$html,&HTML::make_table($col_hdrs,$tab,"Comparing PEGs"));
}

sub show_correlation {
    my($html,$pegs,$vecs) = @_;

    my @stripped = map { ($_ =~ /^fig\|\d+\.\d+\.peg\.(\d+)/) ? $1 : $_ } @$pegs;
    my $col_hdrs = ['PEG',@stripped];
    my $tab = [];
    my $i;
    for ($i=0; ($i < @$pegs); $i++)
    {
	my @corr = map { &corr($vecs,$pegs->[$i],$_) } @$pegs;
	push(@$tab,[&peg_link($pegs->[$i]),@corr]);
    }
    push(@$html,&HTML::make_table($col_hdrs,$tab,'Correlations'));
}


sub show_pearson_correlation {
    my($dir,$html,$pegs) = @_;
    if (-e "$dir/pearson.DB")
    {
	my($pearsonH,$pearson_hash_tieP) = &get_pearson_DB($dir);
	my @stripped = map { ($_ =~ /^fig\|\d+\.\d+\.peg\.(\d+)/) ? $1 : $_ } @$pegs;
	my $col_hdrs = ['PEG',@stripped];
	my $tab = [];
	my $i;
	for ($i=0; ($i < @$pegs); $i++)
	{
	    my @corr = map { &pearson_corr($pearsonH,$pegs->[$i],$_) } @$pegs;
	    push(@$tab,[&peg_link($pegs->[$i]),@corr]);
	}
	push(@$html,&HTML::make_table($col_hdrs,$tab,'Pearson Coefficients'));
	&untie_pearson_hash($pearsonH,$pearson_hash_tieP);
    }
}

sub pearson_corr {
    my($hash,$peg1,$peg2) = @_;
    $peg1 =~ /\.(\d+)$/;
    my $p1 = $1;
    $peg2 =~ /\.(\d+)$/;
    my $p2 = $1;
    my $key = ($p1 < $p2) ? "$p1,$p2" : "$p2,$p1";
    return defined($hash->{$key}) ? $hash->{$key} : " ";
}

sub corr {
    my($vecs,$x,$y) = @_;

    my $v1 = $vecs->{$x};
    my $v2 = $vecs->{$y};
    if (($v1 && $v2) && ($x ne $y))
    {
	my $same = 0;
	my $diff = 0;
	my $i;
	for ($i=0; ($i < @$v1); $i++)
	{
	    if (($v1->[$i] != 1) && ($v2->[$i] ne 1))
	    {
		if ($v1->[$i] == $v2->[$i])
		{
		    $same++;
		}
		else
		{
		    $diff++;
		}
	    }
	}
	if ($same || $diff)
	{
	    return sprintf("%0.2f",(($same - $diff) / ($same + $diff)));
	}
    }
    return " ";
}

sub expressed_val {
    my($v,$val) = @_;

    my $i;
    my $in = [];
    for ($i=0; ($i < @$v); $i++)
    {
	if ($v->[$i] == $val)
	{
	    push(@$in,$i);
	}
    }
    return $in;
}

sub expressed_in {
    my($v) = @_;
    
    return &expressed_val($v,2);
}

sub not_expressed_in {
    my($v) = @_;
    
    return &expressed_val($v,0);
}


MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3