[Bio] / FigWebServices / staph_exp.cgi Repository:
ViewVC logotype

View of /FigWebServices/staph_exp.cgi

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (download) (annotate)
Fri Jun 4 11:35:05 2010 UTC (9 years, 5 months ago) by overbeek
Branch: MAIN
CVS Tags: mgrast_dev_08112011, mgrast_dev_08022011, rast_rel_2014_0912, myrast_rel40, mgrast_dev_05262011, mgrast_dev_04082011, rast_rel_2010_0928, mgrast_version_3_2, mgrast_dev_12152011, mgrast_dev_06072011, rast_rel_2014_0729, 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, mgrast_dev_10262011, HEAD
supports just Staph expression analysis

#########################################################################
# -*- 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 Data::Dumper;
use Statistics::Descriptive;
use tree_utilities;
use HTML;
use strict;

use CGI;
my $cgi = new CGI;

if (0)
{
    my $VAR1;
    eval(join("",`cat $FIG_Config::temp/staph_exp_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/staph_exp_cgi"))
	{
	    print TMP &Dumper($cgi);
	    close(TMP);
	}
    }
    exit;
}
my $html = [];
my $genome        = '158878.1';
my @cel1          = $cgi->param('cel1');
my @cel2          = $cgi->param('cel2');

if (! ((@cel1 > 0)  && (@cel2 > 0)))
{
    &display_initial_page($cgi,$html,$genome);
}
else
{
    &construct_summary($cgi,$html,$genome,\@cel1,\@cel2);
}

&HTML::show_page($cgi,$html);
exit;

sub display_initial_page {
    my($cgi,$html,$genome) = @_;

    push(@$html,$cgi->start_form(-action => "staph_exp.cgi"));
    my @exp = &existing_experiments($genome);
    push(@$html,"<br><br>");
    push(@$html,$cgi->scrolling_list( -name => 'cel1', -values => \@exp, -size => 10, -multiple => 1));
    push(@$html,$cgi->scrolling_list( -name => 'cel2', -values => \@exp, -size => 10, -multiple => 1));
    push(@$html,"<br><br>");
    push(@$html, $cgi->submit('Do It'), $cgi->end_form);
}

sub construct_summary {
    my($cgi,$html,$genome,$cel1,$cel2) = @_;

    my @avg1 = &average_exp($genome,$cel1);  ## returns list of [peg,avg-exp] tuples
    my @avg2 = &average_exp($genome,$cel2);
    my(@diff,$i);
    for ($i=0; ($i < @avg1); $i++)
    {
	($avg1[$i]->[0] eq $avg2[$i]->[0]) || die "incompatible vectors";
	$diff[$i] = [$avg1[$i]->[0],$avg1[$i]->[1] - $avg2[$i]->[1]];
    }
    my @sorted = sort { abs($b->[1]) <=> abs($a->[1]) } @diff;
    my $ids_of_pegs = &affy_ids($genome);
    my $col_hdrs = ['Diff Exp','PEG','Affy','Function','Subsystems'];
    my $fig = new FIG;
    my $tab = [];
    foreach my $tuple (@sorted)
    {
	my($peg,$v) = @$tuple;
	if (abs($v) >= 2)
	{
	    $v = sprintf("%0.3f",$v);
	    my $func    = $fig->function_of($peg);
	    my @subsys = grep {$fig->usable_subsystem($_) } $fig->peg_to_subsystems($peg);
	    push(@$tab,[$v,&peg_link($peg),&affy_id($ids_of_pegs,$peg),$func,join(",",@subsys)]);
	}
    }
    push(@$html,&HTML::make_table($col_hdrs,$tab,"Differential Expression Values"));
}

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

    if ($peg =~ /^fig\|\d+\.\d+\.peg\./)
    {
	my $genome = &SeedUtils::genome_of($peg);
	my $window = "PEG.$$";
	return "<a target=$window href=seedviewer.cgi?genome=$genome&page=Annotation&feature=$peg>$peg</a>";
    }
    else
    {
	return $peg;
    }
}

sub affy_ids {
    my($genome) = @_;

    my $dir = "$FIG_Config::organisms/$genome/UserSpace/ExpressionData";
    my %probe_affy = map { ($_ =~ /^(\d+)\s+(\S+)/) ? ($1 => $2) : () } `cut -f1,2 $dir/probe.info`;

    my %peg_affy;
    foreach $_ (`cut -f1,2 $dir/probepeg.info`)
    {
	if ($_ =~ /^(\d+)\s+(fig\S+)/)
	{
	    $peg_affy{$2}->{$probe_affy{$1}} = 1;
	}
    }
    return \%peg_affy;
}

sub affy_id {
    my($peg_affy,$peg) = @_;
    return join(",",sort keys(%{$peg_affy->{$peg}}));
}

sub average_exp {
    my($genome,$cels) = @_;

    my $v = &raw_vals($genome,$cels->[0]);
    my $i;
    for ($i=1; ($i < @$cels); $i++)
    {
	&add($v,&raw_vals($genome,$cels->[$i]));
    }
    return map { $_->[1] /= @$cels; $_ } @$v;
}

sub raw_vals {
    my($genome,$exp) = @_;

    my $dir = "$FIG_Config::organisms/$genome/UserSpace/ExpressionData";
    my @tmp = `grep \"$exp.CEL.gz\" $dir/experiment.names`;
    my @v;
    if ((@tmp == 1) && ($tmp[0] =~ /^(\d+)/))
    {
	my $col = $1+1;
	@v  = map { chomp; [split(/\s+/,$_)] } `cut -f1,$col $dir/raw_data.tab`;
    }
    else
    {
	@v = ();
    }
    return \@v;
}

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

    my $i;
    for ($i=0; ($i < @$v1); $i++)
    {
	($v1->[$i]->[0] eq $v2->[$i]->[0]) || die "incompatible vectors";
	$v1->[$i]->[1]  += $v2->[$i]->[1];
    }
}

sub existing_experiments {
    my($genome) = @_;

    my $dir           = "$FIG_Config::organisms/$genome/UserSpace/ExpressionData/Experiments";
    opendir(DIR,$dir) || return ();
    my @exp = sort { (lc $a) cmp (lc $b) } map { ($_ =~ /^([^\.].*)\.CEL(\.gz)?$/) ? $1 : () } readdir(DIR);
    closedir(DIR);
    return @exp;
}

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3