[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.2 - (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 : dejongh 1.2 my ($genome_id, $kegg_org_abbrev) = @ARGV;
12 : dejongh 1.1
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 : dejongh 1.2 print ALL_PEGS "peg\tgene_list\tsubsystem\trole\tclassification0\tclassification1\treactions_for_role\tgene_reactions\tkgml_reactions\tmatched_rids2\n";
22 :    
23 : dejongh 1.1 my (%reaction_abbrev_2_ids);
24 :    
25 :     while (my $line = <REACTIONS>)
26 :     {
27 :     chomp($line);
28 :     my ($rabbrev, $rid_list) = split "\t", $line;
29 :     my @rids = split " ", $rid_list;
30 :     $reaction_abbrev_2_ids{$rabbrev} = \@rids;
31 :     }
32 :    
33 :     my (%gene_2_reactions);
34 :    
35 :     while (my $line = <MODEL_GENE_REACTIONS>)
36 :     {
37 :     chomp($line);
38 :     my ($gene, $rabbrev_list) = split "\t", $line;
39 :     my @rabbrevs = split " ", $rabbrev_list;
40 :     foreach my $rabbrev (@rabbrevs)
41 :     {
42 :     if (defined $reaction_abbrev_2_ids{$rabbrev})
43 :     {
44 :     my @rids = @{$reaction_abbrev_2_ids{$rabbrev}};
45 :     push @{$gene_2_reactions{$gene}}, @rids;
46 :     }
47 :     }
48 :     }
49 :    
50 :     my %peg_2_genes;
51 :    
52 :     while (my $line = <MODEL_GENE_PEGS>)
53 :     {
54 :     chomp($line);
55 :     my ($gene, $peg_list) = split "\t", $line;
56 :    
57 :     foreach my $peg (split " ", $peg_list)
58 :     {
59 :     push @{$peg_2_genes{$peg}}, $gene;
60 :     }
61 :     }
62 :    
63 :     foreach my $peg (keys %peg_2_genes)
64 :     {
65 :     my @gene_list = @{$peg_2_genes{$peg}};
66 :     my @gene_reactions = ();
67 : dejongh 1.2 my %kgml_reactions = ();
68 : dejongh 1.1
69 :     foreach my $gene (@gene_list)
70 :     {
71 :     if (defined $gene_2_reactions{$gene})
72 :     {
73 :     push @gene_reactions, @{$gene_2_reactions{$gene}};
74 :     }
75 : dejongh 1.2
76 :     my @tmp = `find $FIG_Config::data/KGML/$kegg_org_abbrev -name "*.xml" -exec grep $gene {} \\;`;
77 :    
78 :     foreach my $line (@tmp)
79 :     {
80 :     if ($line =~ /reaction="(.*)"/)
81 :     {
82 :     my $reactions = $1;
83 :    
84 :     while ($reactions =~ /rn:(R\d+)/g)
85 :     {
86 :     $kgml_reactions{$1} = 1;
87 :     }
88 :     }
89 :     }
90 : dejongh 1.1 }
91 :    
92 :     my $found_ss = 0;
93 :    
94 :     foreach my $ss_info ($fig->subsystems_for_peg($peg))
95 :     {
96 :     my ($subsystem,$role) = @$ss_info;
97 :     my @classification = @{$fig->subsystem_classification($subsystem)};
98 :     my $ss = new Subsystem($subsystem, $fig, 0);
99 :     if (! defined($ss)) { next; }
100 :     $found_ss = 1;
101 :     my $ss_reactions = $ss->get_reactions();
102 :     my $reactions_for_role = $ss_reactions->{$role};
103 :     my %temp = ();
104 : dejongh 1.2 my %intersection1 = ();
105 :     my %intersection2 = ();
106 : dejongh 1.1
107 :     foreach my $rid (@$reactions_for_role) { $temp{$rid} = 1 };
108 :    
109 :     foreach my $rid (@gene_reactions)
110 :     {
111 :     if ($temp{$rid})
112 :     {
113 : dejongh 1.2 $intersection1{$rid} = 1;
114 :     }
115 :     }
116 :    
117 :     my @matched_rids1 = keys %intersection1;
118 :    
119 :     foreach my $rid (keys %kgml_reactions)
120 :     {
121 :     if ($intersection1{$rid})
122 :     {
123 :     $intersection2{$rid} = 1;
124 : dejongh 1.1 }
125 :     }
126 :    
127 : dejongh 1.2 my @matched_rids2 = keys %intersection2;
128 : dejongh 1.1
129 : dejongh 1.2 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@{[ keys %kgml_reactions ]}\t@matched_rids2\n";
130 : dejongh 1.1 }
131 :    
132 :     if (! $found_ss)
133 :     {
134 :     my $func = $fig->function_of($peg);
135 :     print ALL_PEGS "$peg\t@gene_list\t\t$func\t\t\t\t@gene_reactions\t\t\n";
136 :     }
137 :     }
138 :     close(ALL_PEGS);
139 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3