[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.1 - (view) (download) (as text)

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3