Parent Directory
| 
Revision Log
Revision 1.2 - (view) (download) (as text)
| 1 : | formsma | 1.1 | # _*_ Perl _*_ |
| 2 : | # | ||
| 3 : | # | ||
| 4 : | # | ||
| 5 : | use strict; | ||
| 6 : | formsma | 1.2 | use FIGV; |
| 7 : | formsma | 1.1 | use Scenario; |
| 8 : | |||
| 9 : | formsma | 1.2 | my ($genome_id,$orgdir) = @ARGV; |
| 10 : | formsma | 1.1 | |
| 11 : | my $debug = 0; | ||
| 12 : | |||
| 13 : | formsma | 1.2 | |
| 14 : | unless (scalar (@ARGV) >= 1) | ||
| 15 : | formsma | 1.1 | { |
| 16 : | print STDERR "usage: make_SBML_model <genome_id>\n"; | ||
| 17 : | exit(-1); | ||
| 18 : | } | ||
| 19 : | chomp $genome_id; | ||
| 20 : | formsma | 1.2 | my $fig; |
| 21 : | if($orgdir) | ||
| 22 : | { | ||
| 23 : | $fig = new FIGV($orgdir); | ||
| 24 : | } | ||
| 25 : | else | ||
| 26 : | { | ||
| 27 : | $fig = new FIGV; | ||
| 28 : | } | ||
| 29 : | Scenario::set_fig($fig); | ||
| 30 : | formsma | 1.1 | print STDERR "Creating SBML model for genome $genome_id\n" if($debug); |
| 31 : | print STDERR "\tStage 1 - Loading Scenarios, model inputs and biomass\n" if($debug); | ||
| 32 : | my @scenarios = @{Scenario->get_genome_scenarios($genome_id,1)}; | ||
| 33 : | formsma | 1.2 | my $filebase = $fig->model_directory($genome_id)."/Analysis/"; |
| 34 : | formsma | 1.1 | |
| 35 : | my %input_cpds; | ||
| 36 : | open(INPUTS,$filebase."inputs.txt") | ||
| 37 : | or die ("Failed to open ".$filebase."inputs.txt. Run find_model_information $genome_id inputs"); | ||
| 38 : | while(<INPUTS>) | ||
| 39 : | { | ||
| 40 : | chomp; | ||
| 41 : | my @line = split "\t", $_; | ||
| 42 : | chomp @line; | ||
| 43 : | $input_cpds{$line[0]} = 1; | ||
| 44 : | } | ||
| 45 : | close(INPUTS); | ||
| 46 : | #we are going to load up the cofactor compounds into inputs as well | ||
| 47 : | open(COFACT,$filebase."cofactors.txt") | ||
| 48 : | or die("Failed to open ".$filebase."cofactors.txt"); | ||
| 49 : | while(<COFACT>) | ||
| 50 : | { | ||
| 51 : | chomp; | ||
| 52 : | my @line = split "\t",$_; | ||
| 53 : | chomp @line; | ||
| 54 : | $input_cpds{$line[0]} = 1; | ||
| 55 : | } | ||
| 56 : | close(COFACT); | ||
| 57 : | |||
| 58 : | |||
| 59 : | my %compound_id_to_name; | ||
| 60 : | open(IN,"$FIG_Config::global/Models/compound_id_to_name.txt") or die("Failed to open compound_id_to_name.txt"); | ||
| 61 : | while($_ = <IN>){ | ||
| 62 : | chomp($_); | ||
| 63 : | my($id,$name) = split("\t",$_); | ||
| 64 : | $compound_id_to_name{$id} = $name; | ||
| 65 : | } | ||
| 66 : | close(IN); | ||
| 67 : | |||
| 68 : | my %biomass_cpds; | ||
| 69 : | open(BIOMASS,$filebase."biomass_valid.txt") | ||
| 70 : | or die ("Failed to open ".$filebase."biomass_valid.txt\n"); | ||
| 71 : | while(<BIOMASS>) | ||
| 72 : | { | ||
| 73 : | chomp; | ||
| 74 : | my @temp = split "\t", $_; | ||
| 75 : | chomp @temp; | ||
| 76 : | $biomass_cpds{$temp[0]} = 1; | ||
| 77 : | } | ||
| 78 : | close(BIOMASS); | ||
| 79 : | |||
| 80 : | open(IN,$filebase."scenarios_valid.txt"); | ||
| 81 : | my @lines = <IN>; | ||
| 82 : | chomp @lines; | ||
| 83 : | close(IN); | ||
| 84 : | my %valid_scenarios; | ||
| 85 : | foreach (@lines) | ||
| 86 : | { | ||
| 87 : | $valid_scenarios{$_} = 1; | ||
| 88 : | } | ||
| 89 : | |||
| 90 : | my %output_cpds; | ||
| 91 : | open(OUTS,$filebase."outputs.txt"); | ||
| 92 : | while(<OUTS>) | ||
| 93 : | { | ||
| 94 : | chomp; | ||
| 95 : | my @temp = split "\t", $_; | ||
| 96 : | chomp @temp; | ||
| 97 : | $output_cpds{$temp[0]} = 1; | ||
| 98 : | } | ||
| 99 : | close(OUTS); | ||
| 100 : | |||
| 101 : | |||
| 102 : | |||
| 103 : | #load core_metabolites, scenario reaction information | ||
| 104 : | my %core_metabolites; | ||
| 105 : | my %reaction_to_substrates; | ||
| 106 : | my %reaction_to_products; | ||
| 107 : | my %reaction_reversible; | ||
| 108 : | my %reaction_notes; | ||
| 109 : | foreach my $scenario (@scenarios) | ||
| 110 : | { | ||
| 111 : | if(!defined $valid_scenarios{$scenario->get_scenario_name}) | ||
| 112 : | { | ||
| 113 : | next; | ||
| 114 : | } | ||
| 115 : | $scenario->load_reaction_pegs($fig); | ||
| 116 : | my $subsystem = $scenario->get_subsystem_name(); | ||
| 117 : | foreach my $reaction (@{$scenario->get_reaction_ids()}) | ||
| 118 : | { | ||
| 119 : | $reaction_to_substrates{$reaction} = {}; | ||
| 120 : | $reaction_to_products{$reaction} = {}; | ||
| 121 : | |||
| 122 : | my %substrates = %{$scenario->get_reaction_substrates($reaction)}; | ||
| 123 : | my %products = %{$scenario->get_reaction_products($reaction)}; | ||
| 124 : | |||
| 125 : | $reaction_reversible{$reaction} = $scenario->get_reaction_reversibility($reaction); | ||
| 126 : | map {$core_metabolites{$_} = 1; | ||
| 127 : | $reaction_to_substrates{$reaction}->{$_} = $substrates{$_}; | ||
| 128 : | } keys %substrates; | ||
| 129 : | map {$core_metabolites{$_} = 1; | ||
| 130 : | $reaction_to_products{$reaction}->{$_} = $products{$_}; | ||
| 131 : | } keys %products; | ||
| 132 : | if(!defined $reaction_notes{$reaction}) | ||
| 133 : | { | ||
| 134 : | $reaction_notes{$reaction} = {}; | ||
| 135 : | } | ||
| 136 : | foreach ( @{$scenario->get_reaction_pegs($reaction) }) | ||
| 137 : | { | ||
| 138 : | $reaction_notes{$reaction}->{"SEED_PEG:$_"} = 1; | ||
| 139 : | my $roles = $fig->protein_subsystem_to_roles($_, $subsystem); | ||
| 140 : | foreach my $role (@$roles) | ||
| 141 : | { | ||
| 142 : | $reaction_notes{$reaction}->{"Functional Role:$role"} = 1; | ||
| 143 : | } | ||
| 144 : | } | ||
| 145 : | $reaction_notes{$reaction}->{"Subsystem/Scenario:".$scenario->get_scenario_name()} = 1; | ||
| 146 : | |||
| 147 : | } | ||
| 148 : | } | ||
| 149 : | |||
| 150 : | print STDERR "\tStage 1 - Complete\n" if($debug); | ||
| 151 : | |||
| 152 : | print STDERR "\tStage 2 - Write SBML file\n" if($debug); | ||
| 153 : | open(OUT,">".$filebase.$genome_id."_FBA_Model.xml"); | ||
| 154 : | write_header(); | ||
| 155 : | write_species_list(); | ||
| 156 : | write_reactions_list(); | ||
| 157 : | print OUT "</model>\n</sbml>\n"; | ||
| 158 : | close(OUT); | ||
| 159 : | |||
| 160 : | print STDERR "\tStage 2 - Complete\n" if($debug); | ||
| 161 : | |||
| 162 : | |||
| 163 : | sub write_header | ||
| 164 : | { | ||
| 165 : | |||
| 166 : | my $print_id = $genome_id; | ||
| 167 : | $print_id =~ s/\.//g; | ||
| 168 : | print OUT "<?xml version=\"1.0\" encoding=\"UTF-8\"?> | ||
| 169 : | <sbml xmlns=\"http://www.sbml.org/sbml/level2\" level=\"2\" version=\"1\" xmlns:html=\"http://www.w3.org/1999/xhtml\"> | ||
| 170 : | <model id=\"Genome_$print_id\" name=\"Genome_$print_id\"> | ||
| 171 : | <listOfUnitDefinitions> | ||
| 172 : | <unitDefinition id=\"mmol_per_gDW_per_hr\"> | ||
| 173 : | <listOfUnits> | ||
| 174 : | <unit kind=\"mole\" scale=\"-3\"/> | ||
| 175 : | <unit kind=\"gram\" exponent=\"-1\"/> | ||
| 176 : | <unit kind=\"second\" multiplier=\".00027777\" exponent=\"-1\"/> | ||
| 177 : | </listOfUnits> | ||
| 178 : | </unitDefinition> | ||
| 179 : | </listOfUnitDefinitions> | ||
| 180 : | <listOfCompartments> | ||
| 181 : | <compartment id=\"Extra_organism\"/> | ||
| 182 : | <compartment id=\"Cytosol\" outside=\"Extra_organism\"/> | ||
| 183 : | </listOfCompartments>\n"; | ||
| 184 : | } | ||
| 185 : | |||
| 186 : | sub write_species_list | ||
| 187 : | { | ||
| 188 : | print OUT "<listOfSpecies>\n"; | ||
| 189 : | |||
| 190 : | #input compounds | ||
| 191 : | foreach my $id (sort keys %input_cpds) | ||
| 192 : | { | ||
| 193 : | my $name = $compound_id_to_name{$id}; | ||
| 194 : | $name =~ s/\\"//g; | ||
| 195 : | make_species($id,$name,1,1); | ||
| 196 : | } | ||
| 197 : | #biomass compounds | ||
| 198 : | foreach my $id (sort keys(%biomass_cpds)) | ||
| 199 : | { | ||
| 200 : | my $name = $compound_id_to_name{$id}; | ||
| 201 : | $name =~ s/\\"//g; | ||
| 202 : | if(!$input_cpds{$id}) | ||
| 203 : | { | ||
| 204 : | make_species($id,$name,1,1); | ||
| 205 : | } | ||
| 206 : | } | ||
| 207 : | #output compounds | ||
| 208 : | foreach my $id (sort keys(%output_cpds)) | ||
| 209 : | { | ||
| 210 : | my $name = $compound_id_to_name{$id}; | ||
| 211 : | $name =~ s/\\"//g; | ||
| 212 : | if(!$input_cpds{$id} && !$biomass_cpds{$id}) | ||
| 213 : | { | ||
| 214 : | make_species($id,$name,1,1); | ||
| 215 : | } | ||
| 216 : | } | ||
| 217 : | |||
| 218 : | #core_metabolites are all the compounds used in any scenario. | ||
| 219 : | foreach my $id (sort keys(%core_metabolites)) | ||
| 220 : | { | ||
| 221 : | my $name = $compound_id_to_name{$id}; | ||
| 222 : | $name =~ s/\\"//g; | ||
| 223 : | if(!$input_cpds{$id} && !$biomass_cpds{$id} && !$output_cpds{$id}) | ||
| 224 : | { | ||
| 225 : | make_species($id,$name,1,1); | ||
| 226 : | } | ||
| 227 : | } | ||
| 228 : | |||
| 229 : | print OUT "\t<species id=\"biomass_out\" name=\"BIOMASS\" compartment=\"Cytosol\" charge=\"0\" boundaryCondition=\"false\"/>\n"; | ||
| 230 : | |||
| 231 : | print OUT "</listOfSpecies>\n"; | ||
| 232 : | } | ||
| 233 : | |||
| 234 : | sub write_reactions_list | ||
| 235 : | { | ||
| 236 : | print OUT "<listOfReactions>\n"; | ||
| 237 : | |||
| 238 : | #make uptake/transport reactions for the input compounds (provide cpds for the system) | ||
| 239 : | foreach my $input (keys %input_cpds) | ||
| 240 : | { | ||
| 241 : | make_reaction($input."_transport_in",0,{$input."_e" => "1.000000"},{$input => "1.000000"},"0.000000",()); | ||
| 242 : | } | ||
| 243 : | |||
| 244 : | #Make sink/transport reactions for the output compounds | ||
| 245 : | foreach my $output (keys %output_cpds) | ||
| 246 : | { | ||
| 247 : | make_reaction($output."_transport_out",0,{$output => "1.000000"},{$output."_e" => "1.000000"},"0.000000"); | ||
| 248 : | } | ||
| 249 : | |||
| 250 : | #make sinks/transport reactions for all cpds not already outputed or biomassed | ||
| 251 : | foreach my $id (keys %core_metabolites) | ||
| 252 : | { | ||
| 253 : | if(!$output_cpds{$id} && !$biomass_cpds{$id}) | ||
| 254 : | { | ||
| 255 : | make_reaction($id."_sink",0,{$id => "1.000000"},{$id."_e" => "1.000000"},"0.000000"); | ||
| 256 : | } | ||
| 257 : | } | ||
| 258 : | |||
| 259 : | #Now we write out the Scenario reactions | ||
| 260 : | foreach my $rid (keys %reaction_to_substrates) | ||
| 261 : | { | ||
| 262 : | make_reaction($rid,$reaction_reversible{$rid},$reaction_to_substrates{$rid},$reaction_to_products{$rid}, | ||
| 263 : | "0.000000",$reaction_notes{$rid}); | ||
| 264 : | } | ||
| 265 : | |||
| 266 : | #now its time for the biomass reaction | ||
| 267 : | make_reaction("BIOMASS",0,\%biomass_cpds,{"biomass_out" => "1.000000"},"-1.000000"); | ||
| 268 : | |||
| 269 : | print OUT "</listOfReactions>\n"; | ||
| 270 : | } | ||
| 271 : | |||
| 272 : | sub make_species | ||
| 273 : | { | ||
| 274 : | my ($id,$name,$internal,$external) = @_; | ||
| 275 : | if($internal) | ||
| 276 : | { | ||
| 277 : | print OUT "\t<species id=\"$id\" name=\"$name\" compartment=\"Cytosol\" charge=\"0\" boundaryCondition=\"false\"/>\n"; | ||
| 278 : | } | ||
| 279 : | if($external) | ||
| 280 : | { | ||
| 281 : | print OUT "\t<species id=\"$id\_e\" name=\"$name\" compartment=\"Extra_organism\" charge=\"0\" boundaryCondition=\"false\"/>\n"; | ||
| 282 : | } | ||
| 283 : | } | ||
| 284 : | |||
| 285 : | sub make_reaction | ||
| 286 : | { | ||
| 287 : | my($reaction_id,$reversible,$substrates,$products,$objective,$notes) = @_; | ||
| 288 : | |||
| 289 : | |||
| 290 : | if($reversible) | ||
| 291 : | { | ||
| 292 : | print OUT "<reaction id=\"$reaction_id\" reversible=\"true\">\n"; | ||
| 293 : | } | ||
| 294 : | else | ||
| 295 : | { | ||
| 296 : | print OUT "<reaction id=\"$reaction_id\" reversible=\"false\">\n"; | ||
| 297 : | } | ||
| 298 : | if(defined $notes) | ||
| 299 : | { | ||
| 300 : | print OUT "\t<notes>\n"; | ||
| 301 : | foreach my $note (keys %$notes) | ||
| 302 : | { | ||
| 303 : | print OUT "\t\t<html:p>$note</html:p>\n"; | ||
| 304 : | } | ||
| 305 : | print OUT "\t</notes>\n"; | ||
| 306 : | } | ||
| 307 : | print OUT "\t<listOfReactants>\n" if(keys %$substrates); | ||
| 308 : | foreach (keys %$substrates) | ||
| 309 : | { | ||
| 310 : | print OUT "\t\t<speciesReference species=\"$_\" stoichiometry=\"$substrates->{$_}\"/>\n"; | ||
| 311 : | } | ||
| 312 : | print OUT "\t</listOfReactants>\n" if(keys %$substrates); | ||
| 313 : | |||
| 314 : | print OUT "\t<listOfProducts>\n" if(keys %$products); | ||
| 315 : | foreach (keys %$products) | ||
| 316 : | { | ||
| 317 : | print OUT "\t\t<speciesReference species=\"$_\" stoichiometry=\"$products->{$_}\"/>\n"; | ||
| 318 : | } | ||
| 319 : | print OUT "\t</listOfProducts>\n" if(keys %$products); | ||
| 320 : | |||
| 321 : | print OUT "\t<kineticLaw>\n"; | ||
| 322 : | |||
| 323 : | print OUT "\t\t<math xmlns=\"http://www.w3.org/1998/Math/MathML\">\n". | ||
| 324 : | "\t\t\t<apply>\n". | ||
| 325 : | "\t\t\t\t<ci> LOWER_BOUND </ci>\n". | ||
| 326 : | "\t\t\t\t<ci> UPPER_BOUND </ci>\n". | ||
| 327 : | "\t\t\t\t<ci> OBJECTIVE_COEFFICIENT </ci>\n". | ||
| 328 : | "\t\t\t\t<ci> FLUX_VALUE </ci>\n". | ||
| 329 : | "\t\t\t\t<ci> REDUCED_COST </ci>\n". | ||
| 330 : | "\t\t\t</apply>\n\t\t</math>\n"; | ||
| 331 : | |||
| 332 : | |||
| 333 : | print OUT "\t\t<listOfParameters>\n"; | ||
| 334 : | my $low_bound = "0.000000"; | ||
| 335 : | if($reversible) | ||
| 336 : | { | ||
| 337 : | $low_bound = "-999999.000000"; | ||
| 338 : | } | ||
| 339 : | print OUT "\t\t\t<parameter id=\"LOWER_BOUND\" value=\"$low_bound\" units=\"mmol_per_gDW_per_hr\"/>\n" . | ||
| 340 : | "\t\t\t<parameter id=\"UPPER_BOUND\" value=\"999999.000000\" units=\"mmol_per_gDW_per_hr\"/>\n". | ||
| 341 : | "\t\t\t<parameter id=\"OBJECTIVE_COEFFICIENT\" value=\"$objective\"/>\n". | ||
| 342 : | "\t\t\t<parameter id=\"FLUX_VALUE\" value=\"0.000000\" units=\"mmol_per_gDW_per_hr\"/>\n". | ||
| 343 : | "\t\t\t<parameter id=\"REDUCED_COST\" value=\"0.000000\"/>\n"; | ||
| 344 : | print OUT "\t\t</listOfParameters>\n". | ||
| 345 : | "\t</kineticLaw>". | ||
| 346 : | "\n</reaction>\n"; | ||
| 347 : | } |
| MCS Webmaster | ViewVC Help |
| Powered by ViewVC 1.0.3 |