[Bio] / FigKernelScripts / svr_metabolic_reconstruction.pl Repository:
ViewVC logotype

Annotation of /FigKernelScripts/svr_metabolic_reconstruction.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.5 - (view) (download) (as text)

1 : parrello 1.1 #!/usr/bin/perl -w
2 :    
3 : olson 1.4 # This is a SAS Component
4 :    
5 : olson 1.2 use ANNOserver;
6 :     use Getopt::Long;
7 :     use strict;
8 : olson 1.4 use SeedAware;
9 : parrello 1.1
10 :     =head1 svr_metabolic_reconstruction
11 :    
12 :     Get a metabolic reconstruction from a set of functional roles.
13 :    
14 :     By a metabolic reconstruction we mean a detailed list of subsystems
15 :     that we believe are present in a genome or environmental sample, along
16 :     with the specific variant codes that are attached to occurrences of
17 :     the specified subsystems.
18 :    
19 :     If you have a file with tab-delimited fields in which the last field
20 : parrello 1.3 in each line is a functional assignment, then
21 : parrello 1.1
22 :     svr_metabolic_reconstruction < file > extended.file
23 :    
24 :     will produce a file in which each line contains a functional role that
25 : parrello 1.3 is part of an active subsystem. Note that because each functional assignment
26 :     can contain multiple functional roles, a single incoming assignment can produce
27 :     many such lines. For each such line, two extra columns, the variant code and the
28 :     subsystem, will be added. Input lines corresponding to functional roles not
29 :     placed in active subsystems will be written to STDERR (and will not appear in
30 :     STDOUT).
31 : parrello 1.1
32 :     Thus, if one had a FASTA file of protein sequences from some complete
33 :     genome or environmental sample, then
34 :    
35 :     svr_assign_using_figfams < ProteinSequences | svr_metabolic_reconstruction > reconstruction 2> unplaced.hits
36 :    
37 :     would give you the assignments of function plus the unplaced
38 :     sequences.
39 :    
40 :     =cut
41 :    
42 : olson 1.2 my $show_ids;
43 :     my $url = "";
44 :    
45 :     my $usage = "Usage: metabolic_reconstruction [--id] < file";
46 :    
47 :     my $rc = GetOptions("id" => \$show_ids,
48 :     "url=s" => \$url);
49 :    
50 :     if (!$rc || @ARGV > 0)
51 :     {
52 :     die "$usage\n";
53 :     }
54 : parrello 1.1
55 : olson 1.2 my $ss = ANNOserver->new(url => $url);
56 : parrello 1.1
57 : olson 1.5 my $tmp_dir = SeedAware::location_of_tmp();
58 : olson 1.4 my $file = "$tmp_dir/tmpmetabolic.$$";
59 : parrello 1.1
60 :     open(TMP,">$file") || die "could not open $file: $!";
61 :    
62 :     #spool input file to a tmp file. At the end read it back, look up the role in the results hash, write out.
63 :     #if no hit in teh result hash, write to stderr
64 :    
65 : olson 1.2 my %roles;
66 :    
67 :     while (<STDIN>)
68 :     {
69 :     print TMP $_;
70 :     chomp;
71 :     my @in = split "\t";
72 :     my $last = $#in;
73 :     if ($show_ids)
74 :     {
75 :     $roles{$in[$last]} = $in[$last - 1];
76 :     }
77 :     else
78 :     {
79 :     $roles{$in[$last]} = 0;
80 :     }
81 : parrello 1.1 }
82 : olson 1.2 close TMP;
83 : parrello 1.1
84 : olson 1.2 my @role_id;
85 :     while (my($key, $value) = each %roles) {
86 : parrello 1.1 push (@role_id, [$key, $value]);
87 :     }
88 :    
89 : olson 1.2 my @res = $ss->metabolic_reconstruction(-roles => \@role_id);
90 : parrello 1.1 #returns subsys, role, [id]
91 :    
92 : olson 1.2 my %role_hits;
93 : parrello 1.1 foreach my $r (@res) {
94 : olson 1.2 foreach my $res (@$r) {
95 :     my $ss_var = $res->[0];
96 :     my $role = $res->[1];
97 :     $role_hits{$role} = $ss_var;
98 :     if (defined($show_ids)) {
99 :     my $id = $res->[2];
100 :     $roles{$role} = $id;
101 : parrello 1.1 }
102 : olson 1.2 }
103 : parrello 1.1 }
104 :    
105 :     open(TMP,"<$file") || die "could not open $file: $!";
106 :    
107 :     while (<TMP>) {
108 : olson 1.2 chomp;
109 :     my @in = split "\t";
110 :     my $last = $#in;
111 :     my $role = $in[$last];
112 : parrello 1.1 if ($role_hits{$role}) {
113 : olson 1.2 my $ss_var = $role_hits{$role};
114 :     $ss_var =~ /(^.*:)(.+)/;
115 :     my $ss = $1;
116 :     my $var = $2;
117 :     $ss =~ s/:$//;
118 :     print $_;
119 :     print "\t", $ss, "\t", $var, "\n";
120 : parrello 1.1 } else {
121 : olson 1.2 print STDERR $_, "\n";;
122 : parrello 1.1 }
123 :     }
124 : olson 1.2 close(TMP);
125 : parrello 1.1
126 :     unlink $file;
127 :    
128 :     exit;
129 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3