G::SystemsBiology PathwayAlignment
SummaryIncluded librariesPackage variablesSynopsisDescriptionGeneral documentationMethods
Summary
G::Ecell::Pathway - Perl extension for blah blah blah
Package variables
Privates (from "my" definitions)
$path2ec;
$gene2ec;
$orgpath2ec;
$org2ec;
Included modules
G::DB::SDB
G::Messenger
G::Tools::Blast
SelfLoader
SubOpt
Inherit
Exporter
Synopsis
  use G::Ecell::Pathway;
  blah blah blah
Description
Stub documentation for G::Ecell::Pathway was created by h2xs. It looks like the
author of the extension was negligent enough to leave the stub
unedited.
Blah blah blah.
Methods
align_pathway
No description
Code
get_ecvector
No description
Code
get_enzyme_list
No description
Code
load_enzyme2ec
No description
Code
load_gene2ec
No description
Code
Methods description
None available.
Methods code
align_pathwaydescriptionprevnextTop
sub align_pathway {
    require List::Compare;

    opt_default(db=>"org", left_penalty=>0.2, right_penalty=>0.2, output=>"stdout", cutoff=>"10");
    my @argv = opt_get(@_);
    my $query = shift @argv;

    my $db = opt_val("db");
    my $dbref = $org2ec;

    load_enzyme2ec();

    if($db eq 'orgpath'){
	$dbref = $orgpath2ec;
    }elsif($db eq 'path'){
	$dbref = $path2ec;
    }

    my $lp = opt_val("left_penalty");
    my $rp = opt_val("right_penalty");
    my $output = opt_val("output");
    my $cutoff = opt_val("cutoff");

    if($output eq 'stdout'){
	msg_send "\n\nGPACT: G-language Pathway Alignment and Comparison Tool\n";
	msg_send "                   Copyright 2005 G-language Project\n\n";
	msg_send "Database: $db\n";
	msg_send "Left Penalty:  $lp\n";
	msg_send "Right Penalty: $rp\n";
	msg_send "Showing Top $cutoff\n\n";
	msg_send "Query List: ",scalar keys %{$query}, " enzymes\n\n";
	msg_send "Summary:\n-----------------------------------------------------------------\n\n";
    }

    my $result;
    foreach my $entry (keys %{$dbref}){
	my $lc = List::Compare->new([sort keys %{$query}],[sort keys %{$dbref->{$entry}}]);
	my @intersection = $lc->get_intersection();
	my @left = $lc->get_Lonly;
	my @right = $lc->get_Ronly;
	$result->{$entry}->{entry} = $entry;
	$result->{$entry}->{intersection} = scalar @intersection;
	$result->{$entry}->{intersection_a} =\@ intersection;
	$result->{$entry}->{left} = scalar @left;
	$result->{$entry}->{left_a} =\@ left;
	$result->{$entry}->{right} = scalar @right;
	$result->{$entry}->{right_a} =\@ right;
	$result->{$entry}->{identity} = sprintf "%.2f", 
	    scalar(@intersection)/scalar(keys %{$query})*100;
	$result->{$entry}->{score} = 
	    scalar @intersection - (scalar(@left) * $lp + scalar(@right) * $rp); 
    }

    my @order = sort {$result->{$b}->{score} <=> $result->{$a}->{score}} keys %{$result};
    my $return;
    my $i;
    for($i = 0; $i < $cutoff; $i ++){
	$return->[$i] = $result->{$order[$i]};
	last if ($return->[$i]->{intersection} == 0);

	if($output eq 'stdout'){
	    msg_send(sprintf "Hit %3d: %s         score: %3s (%3.2f%s identity)  common: %3d   left: %3d   right:%3d\n",
		     $i + 1, $return->[$i]->{entry}, $return->[$i]->{score},
		     $return->[$i]->{identity}, '%', $return->[$i]->{intersection},
		     $return->[$i]->{left}, $return->[$i]->{right});
	}
    }

    msg_send "\n\nResults:\n-----------------------------------------------------------------\n\n";

    for($i = 0; $i < $cutoff; $i ++){
	$return->[$i] = $result->{$order[$i]};
	last if ($return->[$i]->{intersection} == 0);

	if($output eq 'stdout'){
	    msg_send(sprintf "Hit %3d: %s            score: %s (%.2f%s identity)  left: %d   right:%d\n\n",
		     $i + 1, $return->[$i]->{entry}, $return->[$i]->{score},
		     $return->[$i]->{identity}, '%', $return->[$i]->{left},
		     $return->[$i]->{right});
	    msg_send("    Intersection:\n");
	    my @tmp = @{$return->[$i]->{intersection_a}};
	    while(scalar @tmp){
		msg_send("        ", map {sprintf "%11s", $_} splice(@tmp, 0, 7), "\n");
	    }
	    msg_send("\n    Left only list:\n");
	    my @tmp = @{$return->[$i]->{left_a}};
	    while(scalar @tmp){
		msg_send("        ", map {sprintf "%11s", $_} splice(@tmp, 0, 7), "\n");
	    }
	    msg_send("\n    Right only list:\n");
	    my @tmp = @{$return->[$i]->{right_a}};
	    while(scalar @tmp){
		msg_send("        ", map {sprintf "%11s", $_} splice(@tmp, 0, 8), "\n");
	    }
	    msg_send("\n\n\n");
	}
    }

    if($output eq 'stdout'){
	msg_send("Total $i hits.\n\n");
    }

    return $return;
}
get_ecvectordescriptionprevnextTop
sub get_ecvector {
    my @argv = opt_get(@_);
    my $gb = opt_as_gb(shift @argv);

    load_gene2ec();

    my $flag = 0;
    foreach my $cds ($gb->cds()){
	if(length $gb->{$cds}->{sp_xref}){
	    $flag = 1;
	    last;
	}
    }

    annotate_with_swissprot($gb) unless($flag);

    my $result = {};
    foreach my $cds ($gb->cds()){
	if(length $gb->{$cds}->{sp_xref}){
	    foreach my $ec (keys %{$gene2ec->{$gb->{$cds}->{sp_xref}}}){
		$result->{$ec} ++;
	    }
	}

	if(length $gb->{$cds}->{EC_number}){
	    foreach my $ec (split(/\s+/, $gb->{$cds}->{EC_number})){
		$result->{$ec} ++;
	    }
	}
    }

    return $result;
}
get_enzyme_listdescriptionprevnextTop
sub get_enzyme_list {
    my $name = shift;

    load_enzyme2ec();

    if($name =~ /:/){
	die("ERROR - get_enzyme_list(): organism/pathway combination\" $name\" not available.\n")
	    unless(defined %{$orgpath2ec->{$name}});
	return $orgpath2ec->{$name};
    }elsif($name =~ /^\d+$/){
	die("ERROR - get_enzyme_list(): pathway\" $name\" not available.\n")
	    unless(defined %{$path2ec->{$name}});
	return $path2ec->{$name};
    }else{
	die("ERROR - get_enzyme_list(): organism\" $name\" not available.\n")
	    unless(defined %{$org2ec->{$name}});
	return $org2ec->{$name};
    }
}
load_enzyme2ecdescriptionprevnextTop
sub load_enzyme2ec {
    require List::Compare;
    my $file = shift;

    if(sdb_exists("GLANG_ORGPATH2EC")){
	$orgpath2ec = sdb_load("GLANG_ORGPATH2EC");
	$org2ec = sdb_load("GLANG_ORG2EC");
	$path2ec = sdb_load("GLANG_PATH2EC");
	return 1;
    }

    unless(-e $file){
	system("wget ftp://ftp.genome.ad.jp/pub/kegg/ligand/enzyme -O /tmp/enzyme -q")
	    unless(-e "/tmp/enzyme");
	$file = "/tmp/enzyme";
    }

    my $tmp = 1;
    my $key;
    my $entry;
    open(FILE, $file) || die($!);
    while(<FILE>){
        chomp;

	my $line = '';
	if(/^(\S+?)\s+(.*)/){
	    $key = $1;
	    $line = $2;
	    $entry = $1 if($key eq 'ENTRY' && $line =~ /EC\s+(.*\d)/);
	}elsif(/^\s+(.*)/){
	    $line = $1;
	}

	if($key eq 'PATHWAY' && $line =~ /PATH: .+?(\d+)\s/){
	    $path2ec->{$1}->{$entry} ++;
	}elsif($key eq 'GENES' && $line =~ /(\S+)\:/){
	    $org2ec->{lc($1)}->{$entry} ++;
	}
    }
    close(FILE);

    foreach my $org (keys %{$org2ec}){
	foreach my $path (keys %{$path2ec}){
	    my $lc = List::Compare->new([sort keys %{$org2ec->{$org}}],[sort keys %{$path2ec->{$path}}]);
	    foreach my $enz ($lc->get_intersection()){
		$orgpath2ec->{"$org:$path"}->{$enz} ++;
	    }
	}
    }

    sdb_save($orgpath2ec, "GLANG_ORGPATH2EC");
    sdb_save($org2ec, "GLANG_ORG2EC");
    sdb_save($path2ec, "GLANG_PATH2EC");

    return 1;
}
load_gene2ecdescriptionprevnextTop
sub load_gene2ec {
    my $file = shift;

    if(sdb_exists("GLANG_GENE2EC")){
	$gene2ec = sdb_load("GLANG_GENE2EC");
	return $gene2ec;
    }

    unless(-e $file){
	system("wget http://www.g-language.org/data/gem/gene2ec.txt -O /tmp/gene2ec.txt -q")
	    unless(-e "/tmp/gene2ec.txt");
	$file = "/tmp/gene2ec.txt";
    }

    open(FILE, $file) || die($!);
    while(<FILE>){
        chomp;
        my ($sp, $ec) = split(/\t/, $_, 2);

	foreach my $splitec (split(/\s+/, $ec)){
	    $gene2ec->{$sp}->{$splitec} ++;
	}
    }
    close(FILE);
    
    sdb_save($gene2ec, "GLANG_GENE2EC");
       
    return $gene2ec;
}
General documentation
AUTHORTop
A. U. Thor, a.u.thor@a.galaxy.far.far.away
SEE ALSOTop
perl(1).