[Bio] / DeJonghStuff / map_model_to_subsystems.pl Repository:
ViewVC logotype

Annotation of /DeJonghStuff/map_model_to_subsystems.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : dejongh 1.1 # -*- perl -*-
2 :    
3 :     ###########################################
4 :     use strict;
5 :    
6 :     use FIG;
7 :     use Subsystem;
8 :    
9 :     my $fig = new FIG;
10 :    
11 :     my ($genome_id) = @ARGV;
12 :    
13 :     # map from model's gene-reaction associations to KEGG reaction ids, and check against
14 :     # corresponding functional roles and reactions for genes in subsystems
15 :    
16 :     open(MODEL_GENE_REACTIONS, "<$FIG_Config::global/Models/$genome_id/gene_to_reaction") or die("Run map_model_to_kegg first");
17 :     open(MODEL_GENE_PEGS, "<$FIG_Config::global/Models/$genome_id/gene_to_peg") or die("Run parse_model first");
18 :     open(REACTIONS, "<$FIG_Config::global/Models/$genome_id/reaction_to_kegg_curated");
19 :     open(ALL_PEGS, ">$FIG_Config::global/Models/$genome_id/pegs_to_subsystem");
20 :    
21 :     my (%reaction_abbrev_2_ids);
22 :    
23 :     while (my $line = <REACTIONS>)
24 :     {
25 :     chomp($line);
26 :     my ($rabbrev, $rid_list) = split "\t", $line;
27 :     my @rids = split " ", $rid_list;
28 :     $reaction_abbrev_2_ids{$rabbrev} = \@rids;
29 :     }
30 :    
31 :     my (%gene_2_reactions);
32 :    
33 :     while (my $line = <MODEL_GENE_REACTIONS>)
34 :     {
35 :     chomp($line);
36 :     my ($gene, $rabbrev_list) = split "\t", $line;
37 :     my @rabbrevs = split " ", $rabbrev_list;
38 :     foreach my $rabbrev (@rabbrevs)
39 :     {
40 :     if (defined $reaction_abbrev_2_ids{$rabbrev})
41 :     {
42 :     my @rids = @{$reaction_abbrev_2_ids{$rabbrev}};
43 :     push @{$gene_2_reactions{$gene}}, @rids;
44 :     }
45 :     }
46 :     }
47 :    
48 :     my %peg_2_genes;
49 :    
50 :     while (my $line = <MODEL_GENE_PEGS>)
51 :     {
52 :     chomp($line);
53 :     my ($gene, $peg_list) = split "\t", $line;
54 :    
55 :     foreach my $peg (split " ", $peg_list)
56 :     {
57 :     push @{$peg_2_genes{$peg}}, $gene;
58 :     }
59 :     }
60 :    
61 :     foreach my $peg (keys %peg_2_genes)
62 :     {
63 :     my @gene_list = @{$peg_2_genes{$peg}};
64 :     my @gene_reactions = ();
65 :    
66 :     foreach my $gene (@gene_list)
67 :     {
68 :     if (defined $gene_2_reactions{$gene})
69 :     {
70 :     push @gene_reactions, @{$gene_2_reactions{$gene}};
71 :     }
72 :     }
73 :    
74 :     my $found_ss = 0;
75 :    
76 :     foreach my $ss_info ($fig->subsystems_for_peg($peg))
77 :     {
78 :     my ($subsystem,$role) = @$ss_info;
79 :     my @classification = @{$fig->subsystem_classification($subsystem)};
80 :     my $ss = new Subsystem($subsystem, $fig, 0);
81 :     if (! defined($ss)) { next; }
82 :     $found_ss = 1;
83 :     my $ss_reactions = $ss->get_reactions();
84 :     my $reactions_for_role = $ss_reactions->{$role};
85 :     my %temp = ();
86 :     my %intersection = ();
87 :    
88 :     foreach my $rid (@$reactions_for_role) { $temp{$rid} = 1 };
89 :    
90 :     foreach my $rid (@gene_reactions)
91 :     {
92 :     if ($temp{$rid})
93 :     {
94 :     $intersection{$rid} = 1;
95 :     }
96 :     }
97 :    
98 :     my @matched_rids = keys %intersection;
99 :    
100 :     print ALL_PEGS "$peg\t@gene_list\t$subsystem\t$role\t$classification[0]\t$classification[1]\t@$reactions_for_role\t@gene_reactions\t", scalar @matched_rids, "\t@matched_rids\n";
101 :     }
102 :    
103 :     if (! $found_ss)
104 :     {
105 :     my $func = $fig->function_of($peg);
106 :     print ALL_PEGS "$peg\t@gene_list\t\t$func\t\t\t\t@gene_reactions\t\t\n";
107 :     }
108 :     }
109 :     close(ALL_PEGS);
110 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3