[Bio] / FigKernelPackages / FIGMODELmodel.pm Repository:
ViewVC logotype

Diff of /FigKernelPackages/FIGMODELmodel.pm

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1, Sun Oct 25 21:42:35 2009 UTC revision 1.4, Mon Dec 21 20:05:25 2009 UTC
# Line 13  Line 13 
13          This is the constructor for the FIGMODELmodel object.          This is the constructor for the FIGMODELmodel object.
14  =cut  =cut
15  sub new {  sub new {
16          my ($class,$id,$figmodel) = @_;          my ($class,$figmodel,$id) = @_;
17          my $self->{"_figmodel"}->[0] = $figmodel;          #Error checking first
18          $self->{"_id"}->[0] = $id;          if (!defined($figmodel)) {
19                    print STDERR "FIGMODELmodel->new(undef,".$id."):figmodel must be defined to create a model object!\n";
20                    return undef;
21            }
22            my $self = {_figmodel => $figmodel};
23            bless $self;
24            if (!defined($id)) {
25                    $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,undef):id must be defined to create a model object");
26                    return undef;
27            }
28    
29            #Checking that the id exists
30            my $tbl = $self->figmodel()->database()->GetDBTable("MODELS");
31            if (!defined($tbl)) {
32                    $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,".$id."):could not load MODELS table. Check database!");
33                    return undef;
34            }
35    
36            #If the id is a number, we get the model row by index
37            my $index = $id;
38            if ($id =~ m/^\d+$/) {
39                    $self->{_data} = $tbl->get_row($id);
40            } else {
41                    $self->{_data} = $tbl->get_row_by_key($id,"id");
42                    if (defined($self->{_data})) {
43                            $index = $tbl->row_index($self->{_data});
44                    }
45            }
46            if (!defined($self->{_data})) {
47                    $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,".$id."):could not find specified id in database!");
48                    return undef;
49            }
50            $self->{_index} = $index;
51            $self->figmodel()->{_models}->{$self->id()} = $self;
52            $self->figmodel()->{_models}->{$self->index()} = $self;
53    
54            return $self;
55    }
56    
57    =head3 config
58    Definition:
59            ref::key value = FIGMODELmodel->config(string::key);
60    Description:
61            Trying to avoid using calls that assume configuration data is stored in a particular manner.
62            Call this function to get file paths etc.
63    =cut
64    sub config {
65            my ($self,$key) = @_;
66            return $self->figmodel()->config($key);
67    }
68    
69    =head3 fail
70    Definition:
71            -1 = FIGMODELmodel->fail();
72    Description:
73            Standard return for failed functions.
74    =cut
75    sub fail {
76            my ($self) = @_;
77            return $self->figmodel()->fail();
78    }
79    
80          return bless $self;  =head3 success
81    Definition:
82            1 = FIGMODELmodel->success();
83    Description:
84            Standard return for successful functions.
85    =cut
86    sub success {
87            my ($self) = @_;
88            return $self->figmodel()->success();
89  }  }
90    
91  =head3 figmodel  =head3 figmodel
# Line 28  Line 96 
96  =cut  =cut
97  sub figmodel {  sub figmodel {
98          my ($self) = @_;          my ($self) = @_;
99          return $self->{"_figmodel"}->[0];          return $self->{"_figmodel"};
100    }
101    
102    =head3 fig
103    Definition:
104            FIGMODEL = FIGMODELmodel->fig();
105    Description:
106            Returns a FIG object
107    =cut
108    sub fig {
109            my ($self) = @_;
110    
111            if (!defined($self->{_fig}) && $self->source() !~ /^MGRAST/) {
112                    if ($self->source() =~ /^RAST/) {
113                            $self->{"_fig"} = $self->figmodel()->fig();#$self->genome());
114                    } else {
115                            $self->{"_fig"} = $self->figmodel()->fig();
116                    }
117            }
118    
119            return $self->{"_fig"};
120    }
121    
122    =head3 mgdata
123    Definition:
124            FIGMODEL = FIGMODELmodel->mgdata();
125    Description:
126            Returns a mgrast database object
127    =cut
128    sub mgdata {
129            my ($self) = @_;
130    
131            if (!defined($self->{_mgdata}) && $self->source() =~ /^MGRAST/) {
132                    require MGRAST;
133                    $self->{_mgdata} = $self->figmodel()->mgrast()->Job->get_objects( { 'genome_id' => $self->genome() } )
134            }
135    
136            return $self->{_mgdata};
137  }  }
138    
139  =head3 id  =head3 id
# Line 39  Line 144 
144  =cut  =cut
145  sub id {  sub id {
146          my ($self) = @_;          my ($self) = @_;
147          return $self->{"_id"}->[0];          return $self->{_data}->{id}->[0];
148    }
149    
150    =head3 owner
151    Definition:
152            string = FIGMODELmodel->owner();
153    Description:
154            Returns the username for the model owner
155    =cut
156    sub owner {
157            my ($self) = @_;
158            return $self->{_data}->{owner}->[0];
159    }
160    
161    =head3 index
162    Definition:
163            string = FIGMODELmodel->index();
164    Description:
165            Returns model index
166    =cut
167    sub index {
168            my ($self) = @_;
169            return $self->{_index};
170    }
171    
172    =head3 name
173    Definition:
174            string = FIGMODELmodel->name();
175    Description:
176            Returns the name of the organism or metagenome sample being modeled
177    =cut
178    sub name {
179            my ($self) = @_;
180    
181            if (!defined($self->{_name})) {
182                    my $source = $self->source();
183                    if ($source =~ /^MGRAST/) {
184                            $self->{_name} = "NA";
185                            if (defined($self->mgdata())) {
186                                    $self->{_name} = $self->mgdata()->genome_name;
187                            }
188                    } elsif (defined($self->stats())) {
189                            $self->{_name} = $self->stats()->{'Organism name'}->[0];
190                    } elsif ($source !~ /^RAST/) {
191                            $self->{_name} = $self->fig()->orgname_of_orgid($self->genome());
192                    }
193            }
194    
195            return $self->{_name};
196    }
197    
198    =head3 stats
199    Definition:
200            string = FIGMODELmodel->stats();
201    Description:
202            Returns model stats
203    =cut
204    sub stats {
205            my ($self) = @_;
206    
207            if (!defined($self->{_stats})) {
208                    $self->{_stats} = $self->figmodel()->database()->GetDBTable("MODEL STATS")->get_row_by_key($self->id(),"Model ID");
209            }
210            return $self->{_stats};
211    }
212    
213    =head3 get_reaction_class
214    Definition:
215            string = FIGMODELmodel->get_reaction_class(string::reaction ID);
216    Description:
217            Returns reaction class
218    =cut
219    sub get_reaction_class {
220            my ($self,$reaction,$nohtml) = @_;
221    
222            if (!-e $self->directory()."ReactionClassification-".$self->id().".tbl") {
223                    if (!defined($self->{_reaction_classes})) {
224                            $self->{_reaction_classes} = $self->figmodel()->database()->load_table($self->directory()."ReactionClassification-".$self->id()."-Complete.tbl",";","|",0,["REACTION"]);
225                            if (!defined($self->{_reaction_classes})) {
226                                    return undef;
227                            }
228                    }
229    
230                    my $ClassRow = $self->{_reaction_classes}->get_row_by_key($reaction,"REACTION");
231                    if (defined($ClassRow) && defined($ClassRow->{CLASS})) {
232                            my $class;
233                            if ($ClassRow->{CLASS}->[0] eq "Positive") {
234                                    $class = "Essential =>";
235                            } elsif ($ClassRow->{CLASS}->[0] eq "Negative") {
236                                    $class = "Essential <=";
237                            } elsif ($ClassRow->{CLASS}->[0] eq "Positive variable") {
238                                    $class = "Active =>";
239                            } elsif ($ClassRow->{CLASS}->[0] eq "Negative variable") {
240                                    $class = "Active <=";
241                            } elsif ($ClassRow->{CLASS}->[0] eq "Variable") {
242                                    $class = "Active <=>";
243                            } elsif ($ClassRow->{CLASS}->[0] eq "Blocked") {
244                                    $class = "Inactive";
245                            } elsif ($ClassRow->{CLASS}->[0] eq "Dead") {
246                                    $class = "Disconnected";
247                            }
248                            if (!defined($nohtml) || $nohtml ne "1") {
249                                    $class = "<span title=\"Flux:".$ClassRow->{MIN}->[0]." to ".$ClassRow->{MAX}->[0]."\">".$class."</span>";
250                            }
251                            return $class;
252                    }
253                    return undef;
254            }
255    
256            if (!defined($self->{_reaction_classes})) {
257                    $self->{_reaction_classes} = $self->figmodel()->database()->load_table($self->directory()."ReactionClassification-".$self->id().".tbl",";","|",0,["REACTION"]);
258                    if (!defined($self->{_reaction_classes})) {
259                            return undef;
260                    }
261            }
262    
263            my $ClassRow = $self->{_reaction_classes}->get_row_by_key($reaction,"REACTION");
264            my $classstring = "";
265            if (defined($ClassRow) && defined($ClassRow->{CLASS})) {
266                    for (my $i=0; $i < @{$ClassRow->{CLASS}};$i++) {
267                            if (length($classstring) > 0) {
268                                    $classstring .= "<br>";
269                            }
270                            my $NewClass;
271                            if ($ClassRow->{CLASS}->[$i] eq "Positive") {
272                                    $NewClass = $ClassRow->{MEDIA}->[$i].":Essential =>";
273                                    if (!defined($nohtml) || $nohtml ne "1") {
274                                            $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";
275                                    }
276                            } elsif ($ClassRow->{CLASS}->[$i] eq "Negative") {
277                                    $NewClass = $ClassRow->{MEDIA}->[$i].":Essential <=";
278                                    if (!defined($nohtml) || $nohtml ne "1") {
279                                            $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";
280                                    }
281                            } elsif ($ClassRow->{CLASS}->[$i] eq "Positive variable") {
282                                    $NewClass = $ClassRow->{MEDIA}->[$i].":Active =>";
283                                    if (!defined($nohtml) || $nohtml ne "1") {
284                                            $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";
285                                    }
286                            } elsif ($ClassRow->{CLASS}->[$i] eq "Negative variable") {
287                                    $NewClass = $ClassRow->{MEDIA}->[$i].":Active <=";
288                                    if (!defined($nohtml) || $nohtml ne "1") {
289                                            $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";
290                                    }
291                            } elsif ($ClassRow->{CLASS}->[$i] eq "Variable") {
292                                    $NewClass = $ClassRow->{MEDIA}->[$i].":Active <=>";
293                                    if (!defined($nohtml) || $nohtml ne "1") {
294                                            $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";
295                                    }
296                            } elsif ($ClassRow->{CLASS}->[$i] eq "Blocked") {
297                                    $NewClass = $ClassRow->{MEDIA}->[$i].":Inactive";
298                                    if (!defined($nohtml) || $nohtml ne "1") {
299                                            $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";
300                                    }
301                            } elsif ($ClassRow->{CLASS}->[$i] eq "Dead") {
302                                    $NewClass = $ClassRow->{MEDIA}->[$i].":Disconnected";
303                                    if (!defined($nohtml) || $nohtml ne "1") {
304                                            $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";
305                                    }
306                            }
307                            $classstring .= $NewClass;
308                    }
309            }
310            return $classstring;
311    }
312    
313    =head3 get_biomass
314    Definition:
315            string = FIGMODELmodel->get_biomass();
316    Description:
317            Returns data for the biomass reaction
318    =cut
319    sub get_biomass {
320            my ($self) = @_;
321    
322            if (!defined($self->{_biomass})) {
323                    my $rxntbl = $self->reaction_table();
324                    for (my $i=0; $i < $rxntbl->size(); $i++) {
325                            if ($rxntbl->get_row($i)->{"LOAD"}->[0] =~ m/bio\d\d\d\d\d/) {
326                                    $self->{_biomass} = $rxntbl->get_row($i)->{"LOAD"}->[0];
327                                    last;
328                            }
329                    }
330            }
331    
332            return $self->get_reaction_data($self->{_biomass});
333    }
334    
335    =head3 get_reaction_data
336    Definition:
337            string = FIGMODELmodel->get_reaction_data(string::reaction ID);
338    Description:
339            Returns model reaction data
340    =cut
341    sub get_reaction_data {
342            my ($self,$reaction) = @_;
343    
344            if (!defined($self->reaction_table())) {
345                    return undef;
346            }
347            if ($reaction =~ m/^\d+$/) {
348                    $self->reaction_table()->get_row($reaction);
349            }
350            return $self->reaction_table()->get_row_by_key($reaction,"LOAD");
351    }
352    
353    =head3 get_reaction_id
354    Definition:
355            string = FIGMODELmodel->get_reaction_id(string::reaction ID);
356    Description:
357            Returns model reaction id
358    =cut
359    sub get_reaction_id {
360            my ($self,$reaction) = @_;
361    
362            my $Data = $self->get_reaction_data($reaction);
363            if (!defined($Data) || !defined($Data->{LOAD}->[0])) {
364                    return undef;
365            }
366            return $Data->{LOAD}->[0];
367    }
368    
369    =head3 get_reaction_number
370    Definition:
371            int = FIGMODELmodel->get_reaction_number();
372    Description:
373            Returns the number of reactions in the model
374    =cut
375    sub get_reaction_number {
376            my ($self) = @_;
377    
378            if (!defined($self->reaction_table())) {
379                    return 0;
380            }
381    
382            return $self->reaction_table()->size();
383    }
384    
385    =head3 reaction_table
386    Definition:
387            FIGMODELTable = FIGMODELmodel->reaction_table();
388    Description:
389            Returns FIGMODELTable with the reaction list for the model
390    =cut
391    sub reaction_table {
392            my ($self) = @_;
393    
394            if (!defined($self->{_reaction_data})) {
395                    $self->{_reaction_data} = $self->figmodel()->database()->GetDBModel($self->id());
396            }
397    
398            return $self->{_reaction_data};
399    }
400    
401    =head3 reaction_class_table
402    Definition:
403            FIGMODELTable = FIGMODELmodel->reaction_class_table();
404    Description:
405            Returns FIGMODELTable with the reaction class data, and creates the table file  if it does not exist
406    =cut
407    sub reaction_class_table {
408            my ($self) = @_;
409    
410            if (!defined($self->{_reaction_class_table})) {
411                    if (-e $self->directory()."ReactionClassification-".$self->id().$self->selected_version().".tbl") {
412                            $self->{_reaction_class_table} = $self->figmodel()->database()->load_table($self->directory()."ReactionClassification-".$self->id().$self->selected_version().".tbl",";","|",0,["REACTION","CLASS","MEDIA"]);
413                    } else {
414                            $self->{_reaction_class_table} = FIGMODELTable->new(["REACTION","MEDIA","CLASS","MIN","MAX"],$self->directory()."ReactionClassification-".$self->id().$self->selected_version().".tbl",["REACTION","CLASS","MEDIA"],";","|",undef);
415                    }
416            }
417    
418            return $self->{_reaction_class_table};
419    }
420    
421    =head3 compound_class_table
422    Definition:
423            FIGMODELTable = FIGMODELmodel->compound_class_table();
424    Description:
425            Returns FIGMODELTable with the compound class data, and creates the table file  if it does not exist
426    =cut
427    sub compound_class_table {
428            my ($self) = @_;
429    
430            if (!defined($self->{_compound_class_table})) {
431                    if (-e $self->directory()."CompoundClassification-".$self->id().$self->selected_version().".tbl") {
432                            $self->{_compound_class_table} = $self->figmodel()->database()->load_table($self->directory()."CompoundClassification-".$self->id().$self->selected_version().".tbl",";","|",0,["COMPOUND","CLASS","MEDIA"]);
433                    } else {
434                            $self->{_compound_class_table} = FIGMODELTable->new(["COMPOUND","MEDIA","CLASS","MIN","MAX"],$self->directory()."CompoundClassification-".$self->id().$self->selected_version().".tbl",["COMPOUND","CLASS","MEDIA"],";","|",undef);
435                    }
436            }
437    
438            return $self->{_compound_class_table};
439    }
440    
441    =head3 get_essential_genes
442    Definition:
443            [string::peg ID] = FIGMODELmodel->get_essential_genes(string::media condition);
444    Description:
445            Returns an reference to an array of the predicted essential genes during growth in the input media condition
446    =cut
447    sub get_essential_genes {
448            my ($self,$media) = @_;
449    
450            if (!defined($media)) {
451                    $media = "Complete";
452            }
453            if (!defined($self->{_essential_genes}->{$media})) {
454                    $self->{_essential_genes}->{$media} = undef;
455                    if (-e $self->directory()."EssentialGenes-".$self->id().$self->selected_version()."-".$media.".tbl") {
456                            $self->{_essential_genes}->{$media} = $self->figmodel()->database()->load_single_column_file($self->directory()."EssentialGenes-".$self->id().$self->selected_version()."-".$media.".tbl","");
457                    }
458            }
459    
460            return $self->{_essential_genes}->{$media};
461    }
462    
463    =head3 compound_table
464    Definition:
465            FIGMODELTable = FIGMODELmodel->compound_table();
466    Description:
467            Returns FIGMODELTable with the compound list for the model
468    =cut
469    sub compound_table {
470            my ($self) = @_;
471    
472            if (!defined($self->{_compound_table})) {
473                    $self->{_compound_table} = $self->figmodel()->database()->GetDBModelCompounds($self->id());
474            }
475    
476            return $self->{_compound_table};
477    }
478    
479    =head3 get_compound_data
480    Definition:
481            string = FIGMODELmodel->get_compound_data(string::compound ID);
482    Description:
483            Returns model compound data
484    =cut
485    sub get_compound_data {
486            my ($self,$compound) = @_;
487    
488            if (!defined($self->compound_table())) {
489                    return undef;
490            }
491            if ($compound =~ m/^\d+$/) {
492                    $self->compound_table()->get_row($compound);
493            }
494            return $self->compound_table()->get_row_by_key($compound,"DATABASE");
495    }
496    
497    =head3 genome
498    Definition:
499            string = FIGMODELmodel->genome();
500    Description:
501            Returns model genome
502    =cut
503    sub genome {
504            my ($self) = @_;
505            return $self->{_data}->{genome}->[0];
506    }
507    
508    =head3 source
509    Definition:
510            string = FIGMODELmodel->source();
511    Description:
512            Returns model source
513    =cut
514    sub source {
515            my ($self) = @_;
516            return $self->{_data}->{source}->[0];
517    }
518    
519    =head3 rights
520    Definition:
521            1/0 = FIGMODELmodel->rights(string::username);
522    Description:
523            Returns 1 if the input user can view the model, and zero otherwise
524    =cut
525    sub rights {
526            my ($self,$username) = @_;
527    
528            if ($self->public()) {
529                    return 1;
530            }
531            if (!defined($username) || $username eq "NONE") {
532                    return 0;
533            }
534            if (!defined($self->{_userrights}->{$username})) {
535                    if (defined($self->{_data}->{master})) {
536                            for (my $i=0; $i < @{$self->{_data}->{master}};$i++) {
537                                    $self->{_userrights}->{$self->{_data}->{master}->[$i]} = 1;
538                            }
539                    }
540                    if (defined($self->{_data}->{users})) {
541                            for (my $i=0; $i < @{$self->{_data}->{users}};$i++) {
542                                    $self->{_userrights}->{$self->{_data}->{users}->[$i]} = 1;
543                            }
544                    }
545            }
546            return $self->{_userrights}->{$username};
547    }
548    
549    =head3 public
550    Definition:
551            1/0 = FIGMODELmodel->public();
552    Description:
553            Returns 1 if the model is public, and zero otherwise
554    =cut
555    sub public {
556            my ($self) = @_;
557    
558            if (!defined($self->{_public})) {
559                    $self->{_public} = 0;
560                    if (defined($self->{_data}->{users}->[0]) && $self->{_data}->{users}->[0] eq "all") {
561                            $self->{_public} = 1;
562                    }
563            }
564            return $self->{_public};
565    }
566    
567    =head3 directory
568    Definition:
569            string = FIGMODELmodel->directory();
570    Description:
571            Returns model directory
572    =cut
573    sub directory {
574            my ($self) = @_;
575    
576            if (!defined($self->{_directory})) {
577                    my $userdirectory = "";
578                    if ($self->owner() ne "master") {
579                            $userdirectory = $self->owner()."/";
580                    }
581                    my $source = $self->source();
582                    if ($source =~ /^MGRAST/) {
583                            $self->{_directory} = $self->figmodel()->config("organism directory")->[0].$userdirectory.$self->genome()."/";
584                    } elsif ($source =~ /^RAST/) {
585                            $self->{_directory} = $self->figmodel()->config("organism directory")->[0].$userdirectory.$self->genome()."/";
586                    } elsif ($source =~ /^SEED/) {
587                            $self->{_directory} = $self->figmodel()->config("organism directory")->[0].$userdirectory.$self->genome()."/";
588                    } elsif ($source =~ /^PM/) {
589                            if (length($userdirectory) == 0) {
590                                    $self->{_directory} = $self->figmodel()->config("imported model directory")->[0].$self->id()."/";
591                            } else {
592                                    $self->{_directory} = $self->figmodel()->config("organism directory")->[0].$userdirectory.$self->id()."/";
593                            }
594                    }
595            }
596    
597            return $self->{_directory};
598    }
599    
600    =head3 filename
601    Definition:
602            string = FIGMODELmodel->filename();
603    Description:
604            Returns model filename
605    =cut
606    sub filename {
607            my ($self) = @_;
608    
609            return $self->directory().$self->id().$self->selected_version().".txt";
610    }
611    
612    =head3 set_metagenome_stats
613    Definition:
614            string = FIGMODELmodel->set_metagenome_stats();
615    Description:
616            Sets the values of many model stats for a metagenome
617    =cut
618    sub set_metagenome_stats {
619            my ($self) = @_;
620    
621            $self->{_total_compounds} = 0;
622            if (defined($self->compound_table())) {
623                    $self->{_total_compounds} = $self->compound_table()->size();
624            }
625            $self->{_gene_reactions} = 0;
626            $self->{_gapfilling_reactions} = 0;
627            $self->{_model_genes} = 0;
628            $self->{_total_reactions} = 0;
629            if (defined($self->reaction_table())) {
630                    $self->{_total_reactions} = $self->reaction_table()->size();
631                    my $tbl = $self->reaction_table();
632                    my $spontaneous = 0;
633                    for (my $i=0; $i < $tbl->size(); $i++) {
634                            my $row = $tbl->get_row($i);
635                            if (!defined($row->{"ASSOCIATED PEG"}->[0]) || $row->{"ASSOCIATED PEG"}->[0] !~ m/peg/) {
636                                    if ($row->{"ASSOCIATED PEG"}->[0] =~ m/SPONTANEOUS/) {
637                                            $spontaneous++;
638                                    } else {
639                                            $self->{_gapfilling_reactions}++;
640                                    }
641                            } else {
642                                    for (my $j=0; $j < @{$row->{"CONFIDENCE"}}; $j++) {
643                                            my @ecores = split(/;/,$row->{"CONFIDENCE"}->[$j]);
644                                            $self->{_model_genes} += @ecores;
645                                    }
646                            }
647                    }
648                    $self->{_gene_reactions} = $tbl->size() - $spontaneous - $self->{_gapfilling_reactions};
649            }
650    }
651    
652    =head3 version
653    Definition:
654            string = FIGMODELmodel->version();
655    Description:
656            Returns the version of the model
657    =cut
658    sub version {
659            my ($self) = @_;
660    
661            if (!defined($self->{_version})) {
662                    if (!defined($self->{_selectedversion})) {
663                            if (defined($self->stats())) {
664                                    $self->{_version} = "V".$self->stats()->{"Version"}->[0].".".$self->stats()->{"Gap fill version"}->[0];
665                            }
666                    } else {
667                            $self->{_version} = $self->{_selectedversion};
668                    }
669            }
670            return $self->{_version};
671    }
672    
673    =head3 selected_version
674    Definition:
675            string = FIGMODELmodel->selected_version();
676    Description:
677            Returns the selected version of the model
678    =cut
679    sub selected_version {
680            my ($self) = @_;
681    
682            if (!defined($self->{_selectedversion})) {
683                    return "";
684            }
685            return $self->{_selectedversion};
686    }
687    
688    =head3 modification_time
689    Definition:
690            string = FIGMODELmodel->modification_time();
691    Description:
692            Returns the selected version of the model
693    =cut
694    sub modification_time {
695            my ($self) = @_;
696            if (!defined($self->{_modification_time})) {
697                    my $stats = $self->stats();
698                    if (defined($stats)) {
699                            $self->{_modification_time} = 0;
700                            if (defined($stats->{"Build date"}->[0]) && $self->{_modification_time} < $stats->{"Build date"}->[0]) {
701                                    $self->{_modification_time} = $stats->{"Build date"}->[0];
702                            } elsif (defined($stats->{"Gap fill date"}->[0]) && $self->{_modification_time} < $stats->{"Gap fill date"}->[0]) {
703                                    $self->{_modification_time} = $stats->{"Gap fill date"}->[0];
704                            }
705                    } else {
706                            $self->{_modification_time} = 0;
707                    }
708            }
709            return $self->{_modification_time};
710    }
711    
712    =head3 gene_reactions
713    Definition:
714            string = FIGMODELmodel->gene_reactions();
715    Description:
716            Returns the number of reactions added by the gap filling
717    =cut
718    sub gene_reactions {
719            my ($self) = @_;
720    
721            if (!defined($self->{_gene_reactions})) {
722                    if ($self->source() =~ /^MGRAST/) {
723                            $self->set_metagenome_stats();
724                    } elsif (defined($self->stats())) {
725                            $self->{_gene_reactions} = $self->total_reactions() - $self->gapfilling_reactions() - $self->stats()->{'Spontaneous'}->[0] - $self->stats()->{'Growmatch reactions'}->[0] - $self->stats()->{'Biolog gap filling reactions'}->[0];
726                    }
727            }
728            return $self->{_gene_reactions};
729    }
730    
731    =head3 total_compounds
732    Definition:
733            string = FIGMODELmodel->total_compounds();
734    Description:
735            Returns the number of compounds in the model
736    =cut
737    sub total_compounds {
738            my ($self) = @_;
739    
740            if (!defined($self->{_total_compounds})) {
741                    if ($self->source() =~ /^MGRAST/) {
742                            $self->set_metagenome_stats();
743                    } elsif (defined($self->stats())) {
744                            $self->{_total_compounds} = $self->stats()->{'Metabolites'}->[0];
745                    }
746            }
747            return $self->{_total_compounds};
748    }
749    
750    =head3 gapfilling_reactions
751    Definition:
752            string = FIGMODELmodel->gapfilling_reactions();
753    Description:
754            Returns the number of reactions added by the gap filling
755    =cut
756    sub gapfilling_reactions {
757            my ($self) = @_;
758    
759            if (!defined($self->{_gapfilling_reactions})) {
760                    if ($self->source() =~ /^MGRAST/) {
761                            $self->set_metagenome_stats();
762                    } elsif (defined($self->stats())) {
763                            $self->{_gapfilling_reactions} = $self->stats()->{'Gap filling reactions'}->[0];
764                    }
765            }
766            return $self->{_gapfilling_reactions};
767    }
768    
769    =head3 total_reactions
770    Definition:
771            string = FIGMODELmodel->total_reactions();
772    Description:
773            Returns the total number of reactions in the model
774    =cut
775    sub total_reactions {
776            my ($self) = @_;
777    
778            if (!defined($self->{_total_reactions})) {
779                    if ($self->source() =~ /^MGRAST/) {
780                            $self->set_metagenome_stats();
781                    } elsif (defined($self->stats())) {
782                            $self->{_total_reactions} = $self->stats()->{'Number of reactions'}->[0];
783                    }
784            }
785            return $self->{_total_reactions};
786    }
787    
788    =head3 model_genes
789    Definition:
790            string = FIGMODELmodel->model_genes();
791    Description:
792            Returns the number of genes mapped to one or more reactions in the model
793    =cut
794    sub model_genes {
795            my ($self) = @_;
796    
797            if (!defined($self->{_model_genes})) {
798                    if ($self->source() =~ /^MGRAST/) {
799                            $self->set_metagenome_stats();
800                    } elsif (defined($self->stats())) {
801                            $self->{_model_genes} = $self->stats()->{'Genes with reactions'}->[0];
802                    }
803            }
804            return $self->{_model_genes};
805    }
806    
807    =head3 class
808    Definition:
809            string = FIGMODELmodel->class();
810    Description:
811            Returns the class of the model: gram positive, gram negative, other
812    =cut
813    sub class {
814            my ($self) = @_;
815    
816            if (!defined($self->{_class})) {
817                    if ($self->source() =~ /^MGRAST/) {
818                            $self->{_class} = "Other";
819                    } elsif (defined($self->stats())) {
820                            $self->{_class} = $self->stats()->{Class}->[0];
821                    }
822            }
823            return $self->{_class};
824    }
825    
826    =head3 taxonomy
827    Definition:
828            string = FIGMODELmodel->taxonomy();
829    Description:
830            Returns model taxonomy or biome if this is an metagenome model
831    =cut
832    sub taxonomy {
833            my ($self) = @_;
834    
835            if (!defined($self->{_taxonomy})) {
836                    my $source = $self->source();
837                    if ($source =~ /^MGRAST/) {
838                            $self->{_taxonomy} = "NA";
839                            my $mgmodeldata = $self->figmodel()->database()->mg_model_data($self->genome());
840                            if (defined($mgmodeldata)) {
841                                    $self->{_taxonomy} = $mgmodeldata->biome();
842                            }
843                    } else {
844                            $self->{_taxonomy} = $self->fig()->taxonomy_of($self->genome());
845                    }
846            }
847    
848            return $self->{_taxonomy};
849    }
850    
851    =head3 genome_size
852    Definition:
853            string = FIGMODELmodel->genome_size();
854    Description:
855            Returns size of the modeled genome in KB
856    =cut
857    sub genome_size {
858            my ($self) = @_;
859    
860            if (!defined($self->{_genome_size})) {
861                    my $source = $self->source();
862                    if ($source =~ /^MGRAST/) {
863                            $self->{_genome_size} = "NA";
864                            if (defined($self->mgdata())) {
865                                    $self->{_genome_size} = $self->mgdata()->size;
866                            }
867                    } else {
868                            $self->{_genome_size} = $self->fig()->genome_szdna($self->genome());
869                    }
870            }
871    
872            return $self->{_genome_size};
873    }
874    
875    =head3 genome_genes
876    Definition:
877            string = FIGMODELmodel->genome_genes();
878    Description:
879            Returns the number of genes in the modeled genome
880    =cut
881    sub genome_genes {
882            my ($self) = @_;
883    
884            if (!defined($self->{_genome_genes})) {
885                    my $source = $self->source();
886                    if ($source =~ /^MGRAST/) {
887                            $self->{_genome_genes} = "NA";
888                            if (defined($self->mgdata())) {
889                                    $self->{_genome_genes} = $self->mgdata()->genome_contig_count;
890                            }
891                    } elsif (defined($self->stats())) {
892                            $self->{_genome_genes} = $self->stats()->{'Total genes'}->[0];
893                    }
894            }
895    
896            return $self->{_genome_genes};
897    }
898    
899    =head3 run_default_model_predictions
900    Definition:
901            0/1::status = FIGMODELmodel->run_default_model_predictions(string::media ID);
902    Description:
903    =cut
904    sub run_default_model_predictions {
905            my ($self,$Media) = @_;
906    
907            #Assuming complete media if none is provided
908            if (!defined($Media)) {
909                    $Media = "Complete";
910            }
911    
912            #Predicting essentiality
913            my $result = $self->figmodel()->RunFBASimulation($self->id(),"SINGLEKO",undef,undef,[$self->id()],[$Media]);
914            #Checking that the table is defined and the output file exists
915            if (defined($result) && defined($result->get_row(0)->{"ESSENTIALGENES"})) {
916                    $self->figmodel()->database()->print_array_to_file($self->directory()."EssentialGenes-".$self->id()."-".$Media.".tbl",[join("\n",@{$result->get_row(0)->{"ESSENTIALGENES"}})]);
917            } else {
918                    $self->figmodel()->error_message("FIGMODELmodel:run_default_model_predictions:could not identify essential reactions for model ".$self->id().$self->selected_version().".");
919                    return $self->figmodel()->fail();
920            }
921    
922            #Classifying reactions and compounds
923            my $tbl = $self->classify_model_reactions($Media);
924            if (!defined($tbl)) {
925                    $self->figmodel()->error_message("FIGMODELmodel:run_default_model_predictions:could not classify reactions for model ".$self->id().$self->selected_version().".");
926                    return $self->figmodel()->fail();
927            }
928            $tbl->save();
929    
930            return $self->figmodel()->success();
931    }
932    
933    =head3 save_obsolete_stats
934    Definition:
935            FIGMODELmodel->save_obsolete_stats();
936    Description:
937    =cut
938    sub save_obsolete_stats {
939            my ($self) = @_;
940    
941            #checking if stats exists
942            my $stats = $self->stats();
943            if (defined($stats)) {
944                    $stats->{"Model ID"}->[0] = $self->id()."V".$stats->{"Version"}->[0].".".$stats->{"Gap fill version"}->[0];
945                    $self->figmodel()->database()->update_row("OBSOLETE MODEL STATS",$stats,"Model ID");
946                    $stats->{"Model ID"}->[0] = $self->id();
947            }
948    }
949    
950    =head3 update_stats_for_gap_filling
951    Definition:
952            {string => [string]} = FIGMODELmodel->update_stats_for_gap_filling(int::gapfill time);
953    Description:
954    =cut
955    sub update_stats_for_gap_filling {
956            my ($self,$gapfilltime) = @_;
957    
958            #preserving the stats for the now obselete model
959            $self->save_obsolete_stats();
960            my $stats = $self->update_model_stats(0);
961            $stats->{"Gap filling time"}->[0] = $gapfilltime;
962            $stats->{"Gap fill date"}->[0] = time();
963            if (!defined($stats->{"Gap fill version"}->[0])) {
964                    $stats->{"Gap fill version"}->[0] = 0;
965            }
966            $stats->{"Gap fill version"}->[0]++;
967    
968            #Updating the stats stored in the table
969            $self->figmodel()->database()->update_row("MODEL STATS",$stats,"Model ID");
970            return $stats;
971    }
972    
973    =head3 update_stats_for_build
974    Definition:
975            {string => [string]} = FIGMODELmodel->update_stats_for_build();
976    Description:
977    =cut
978    sub update_stats_for_build {
979            my ($self) = @_;
980    
981            #preserving the stats for the now obselete model
982            $self->save_obsolete_stats();
983            my $stats = $self->update_model_stats(0);
984            $stats->{"Build date"}->[0] = time();
985            if (!defined($stats->{"Version"}->[0])) {
986                    $stats->{"Version"}->[0] = 0;
987            }
988            $stats->{"Version"}->[0]++;
989    
990            #Updating the stats stored in the table
991            $self->figmodel()->database()->update_row("MODEL STATS",$stats,"Model ID");
992            return $stats;
993    }
994    
995    =head3 update_model_stats
996    Definition:
997            FIGMODELmodel->update_model_stats();
998    Description:
999    =cut
1000    sub update_model_stats {
1001            my ($self) = @_;
1002    
1003            #Getting reaction table
1004            my $rxntbl = $self->reaction_table();
1005            if (!defined($rxntbl)) {
1006                    $self->figmodel()->error_message("FIGMODELmodel:update_model_stats:Could not load reaction list for ".$self->id());
1007                    return undef;
1008            }
1009            my $cpdtbl = $self->compound_table();
1010    
1011            #Creating empty status row
1012            my $CurrentStats = {"Genes with reactions" => [0],
1013                                                     "Metabolites" => [$cpdtbl->size()],
1014                                                     "Growmatch reactions" => [0],
1015                                                     "Spontaneous" => [0],
1016                                                     "Biolog gap filling reactions" => [0],
1017                                                     "Number of reactions" => [$rxntbl->size()],
1018                                                     "Transport reaction"=>[0],
1019                                                     "Gap filling reactions" => [0],
1020                                                     "Model ID" => [$self->id()],
1021                                                     "Subsystem genes with reactions" => [0],
1022                                                     "Nonsubsystem genes with reactions" => [0],
1023                                                     Version => [0],
1024                                                     "Gap fill version" => [0],
1025                                                     Source => [$self->source()],
1026                                                     "Genes with one reaction" => [0],
1027                                                     "Genome ID" => [$self->genome()]};
1028    
1029            my $genomestats = $self->figmodel()->get_genome_stats($self->genome());
1030            if (defined($genomestats)) {
1031                    $CurrentStats->{SOURCE}->[0] = $genomestats->{SOURCE}->[0];
1032                    $CurrentStats->{"Organism name"}->[0] = $genomestats->{NAME}->[0];
1033                    $CurrentStats->{"Total genes"}->[0] = $genomestats->{"TOTAL GENES"}->[0];
1034                    $CurrentStats->{"Gram positive genes"}->[0] = $genomestats->{"GRAM POSITIVE GENES"}->[0];
1035                    $CurrentStats->{"Gram negative genes"}->[0] = $genomestats->{"GRAM NEGATIVE GENES"}->[0];
1036                    $CurrentStats->{"Class"}->[0] = $genomestats->{CLASS}->[0];
1037                    $CurrentStats->{"Genes with functions"}->[0] = $genomestats->{"GENES WITH FUNCTIONS"}->[0];
1038                    $CurrentStats->{"Subsystem genes"}->[0] = $genomestats->{"SUBSYSTEM GENES"}->[0];
1039                    $CurrentStats->{"Nonsubsystem genes"}->[0] = $genomestats->{"NON SUBSYSTEM GENES"}->[0];
1040                    if (defined($genomestats->{TAXONOMY})) {
1041                            $CurrentStats->{"Taxonomy 0"}->[0] = $genomestats->{TAXONOMY}->[0];
1042                            $CurrentStats->{"Taxonomy 1"}->[0] = $genomestats->{TAXONOMY}->[1];
1043                            $CurrentStats->{"Taxonomy 2"}->[0] = $genomestats->{TAXONOMY}->[2];
1044                            $CurrentStats->{"Taxonomy 3"}->[0] = $genomestats->{TAXONOMY}->[3];
1045                            $CurrentStats->{"Taxonomy 4"}->[0] = $genomestats->{TAXONOMY}->[4];
1046                            $CurrentStats->{"Taxonomy 5"}->[0] = $genomestats->{TAXONOMY}->[5];
1047                    }
1048            }
1049    
1050            #Transfering build, version, and gap fill data from existing stats
1051            if (defined($self->stats())) {
1052                    $CurrentStats->{Version}->[0] = $self->stats()->{Version}->[0];
1053                    $CurrentStats->{"Gap fill version"}->[0] = $self->stats()->{"Gap fill version"}->[0];
1054                    $CurrentStats->{"Gap filling time"}->[0] = $self->stats()->{"Gap filling time"}->[0];
1055                    $CurrentStats->{"Gap fill date"}->[0] = $self->stats()->{"Gap fill date"}->[0];
1056                    $CurrentStats->{"Build date"}->[0] = $self->stats()->{"Build date"}->[0];
1057            }
1058    
1059            my %GeneHash;
1060            my %NonpegHash;
1061            my %CompoundHash;
1062            for (my $i=0; $i < $rxntbl->size(); $i++) {
1063                    my $Row = $rxntbl->get_row($i);
1064                    if (defined($Row) && defined($Row->{"ASSOCIATED PEG"})) {
1065                            my $ReactionRow = $self->figmodel()->get_reaction($Row->{"LOAD"}->[0]);
1066                            if (defined($ReactionRow->{"EQUATION"}->[0])) {
1067                                    #Checking for extracellular metabolites which indicate that this is a transporter
1068                                    if ($ReactionRow->{"EQUATION"}->[0] =~ m/\[e\]/) {
1069                                            $CurrentStats->{"Transport reaction"}->[0]++;
1070                                    }
1071                            }
1072                            #Identifying spontaneous/biolog/gapfilling/gene associated reactions
1073                            if ($Row->{"ASSOCIATED PEG"}->[0] =~ m/BIOLOG/i) {
1074                                    $CurrentStats->{"Biolog gap filling reactions"}->[0]++;
1075                            } elsif ($Row->{"ASSOCIATED PEG"}->[0] =~ m/SPONTANEOUS/i) {
1076                                    $CurrentStats->{"Spontaneous"}->[0]++;
1077                            } elsif ($Row->{"ASSOCIATED PEG"}->[0] =~ m/GAP/ || $Row->{"ASSOCIATED PEG"}->[0] =~ m/UNIVERSAL/i || $Row->{"ASSOCIATED PEG"}->[0] =~ m/UNKNOWN/i) {
1078                                    $CurrentStats->{"Gap filling reactions"}->[0]++;
1079                            } else {
1080                                    foreach my $GeneSet (@{$Row->{"ASSOCIATED PEG"}}) {
1081                                            my @GeneList = split(/\+/,$GeneSet);
1082                                            foreach my $Gene (@GeneList) {
1083                                                    if ($Gene =~ m/(peg\.\d+)/) {
1084                                                            $GeneHash{$1} = 1;
1085                                                    } else {
1086                                                            $NonpegHash{$Gene} = 1;
1087                                                    }
1088                                            }
1089                                    }
1090                            }
1091                    }
1092            }
1093            my @genes = keys(%GeneHash);
1094            my @othergenes = keys(%NonpegHash);
1095            $CurrentStats->{"Genes with reactions"}->[0] = @genes + @othergenes;
1096    
1097            #Updating the stats stored in the table
1098            $self->figmodel()->database()->update_row("MODEL STATS",$CurrentStats,"Model ID");
1099            $self->{_stats} = $CurrentStats;
1100            return $CurrentStats;
1101    }
1102    
1103    =head3 GapFillModel
1104    Definition:
1105            (success/fail) = FIGMODELmodel->GapFillModel();
1106    Description:
1107            This function performs an optimization to identify the minimal set of reactions that must be added to a model in order for biomass to be produced by the biomass reaction in the model.
1108            Before running the gap filling, the existing model is backup in the same directory with the current version numbers appended.
1109            If the model has been gap filled previously, the previous gap filling reactions are removed prior to running the gap filling again.
1110    =cut
1111    sub GapFillModel {
1112            my ($self,$donotclear) = @_;
1113    
1114            #Setting status of model to gap filling
1115            my $OrganismID = $self->genome();
1116            $self->set_status(1,"Auto completion running");
1117            my $UniqueFilename = $self->figmodel()->filename();
1118            my $StartTime = time();
1119    
1120            #Archiving the existing model
1121            $self->ArchiveModel();
1122    
1123            #Removing any gapfilling reactions that may be currently present in the model
1124            if (!defined($donotclear) || $donotclear != 1) {
1125                    my $ModelTable = $self->reaction_table();
1126                    for (my $i=0; $i < $ModelTable->size(); $i++) {
1127                            if (!defined($ModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0]) || $ModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0] eq "GAP FILLING" || $ModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0] =~ m/BIOLOG/ || $ModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0] =~ m/GROWMATCH/) {
1128                                    $ModelTable->delete_row($i);
1129                                    $i--;
1130                            }
1131                    }
1132                    $ModelTable->save();
1133            }
1134    
1135            #Calling the MFAToolkit to run the gap filling optimization
1136            my $MinimalMediaTable = $self->figmodel()->database()->GetDBTable("MINIMAL MEDIA TABLE");
1137            my $Media = "Complete";
1138            if (defined($MinimalMediaTable->get_row_by_key($self->genome(),"Organism"))) {
1139                    $Media = $MinimalMediaTable->get_row_by_key($self->genome(),"Organism")->{"Minimal media"}->[0];
1140                    #Loading media, changing bounds, saving media as a test media
1141                    my $MediaTable = FIGMODELTable::load_table($self->config("Media directory")->[0].$Media.".txt",";","",0,["VarName"]);
1142                    for (my $i=0; $i < $MediaTable->size(); $i++) {
1143                            if ($MediaTable->get_row($i)->{"Min"}->[0] < 0) {
1144                                    $MediaTable->get_row($i)->{"Min"}->[0] = -10000;
1145                            }
1146                            if ($MediaTable->get_row($i)->{"Max"}->[0] > 0) {
1147                                    $MediaTable->get_row($i)->{"Max"}->[0] = 10000;
1148                            }
1149                    }
1150                    $MediaTable->save($self->config("Media directory")->[0].$UniqueFilename."TestMedia.txt");
1151                    system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$UniqueFilename."TestMedia",["GapFilling"],{"Default max drain flux" => 0},"GapFill".$self->id().".log",undef));
1152                    unlink($self->config("Media directory")->[0].$UniqueFilename."TestMedia.txt");
1153            } else {
1154                    system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),undef,["GapFilling"],undef,"GapFill".$self->id().".log",undef));
1155            }
1156    
1157            #Parse the solutions file for the model and get the reaction list from it
1158            my $SolutionData = $self->figmodel()->LoadProblemReport($UniqueFilename);
1159            #Clearing the mfatoolkit output and log file
1160            $self->figmodel()->clearing_output($UniqueFilename,"GapFill".$self->id().".log");
1161            if (!defined($SolutionData) || $SolutionData->size() == 0) {
1162                    $self->set_status(1,"Auto completion failed: auto completion solution not found");
1163                    $self->figmodel()->error_message("FIGMODEL:GapFillModel: Gap filling report not found!");
1164                    return $self->fail();
1165            }
1166            $SolutionData->save("/home/chenry/Solution.txt");
1167    
1168            #Looking for the last printed recursive MILP solution
1169            for (my $i=($SolutionData->size()-1); $i >=0; $i--) {
1170                    if (defined($SolutionData->get_row($i)->{"Notes"}) && $SolutionData->get_row($i)->{"Notes"}->[0] =~ m/^Recursive/) {
1171                            my $AllSolutions = substr($SolutionData->get_row($i)->{"Notes"}->[0],15);
1172                            my @TempThree = split(/\|/,$AllSolutions);
1173                            if (@TempThree > 0 && $TempThree[0] =~ m/.+:(.+)/) {
1174                                    my @TempFour = split(/,/,$1);
1175                                    my $DirectionList;
1176                                    my $ReactionList;
1177                                    for (my $j=0; $j < @TempFour; $j++) {
1178                                            my $ID = "";
1179                                            my $Sign = "=>";
1180                                            if ($TempFour[$j] =~ m/\-(rxn\d\d\d\d\d)/) {
1181                                                    $ID = $1;
1182                                                    $Sign = "<=";
1183                                            } elsif ($TempFour[$j] =~ m/\+(rxn\d\d\d\d\d)/) {
1184                                                    $ID = $1;
1185                                            }
1186                                            if ($ID ne "") {
1187                                                    if ($self->figmodel()->reversibility_of_reaction($ID) ne $Sign) {
1188                                                            $Sign = "<=>";
1189                                                    }
1190                                                    push(@{$DirectionList},$Sign);
1191                                                    push(@{$ReactionList},$ID);
1192                                            }
1193                                    }
1194                                    $self->figmodel()->IntegrateGrowMatchSolution($self->id(),undef,$ReactionList,$DirectionList,"GAP FILLING",0,1);
1195                            }
1196                            last;
1197                    }
1198            }
1199            #Updating model stats with gap filling results
1200            my $ElapsedTime = time() - $StartTime;
1201            $self->figmodel()->ClearDBModel($self->id(),1);
1202            #Determining why each gap filling reaction was added
1203            $self->figmodel()->IdentifyDependancyOfGapFillingReactions($self->id(),$Media);
1204            if (!defined($donotclear) || $donotclear != 1) {
1205                    $self->figmodel()->AddBiologTransporters($self->id());
1206                    if ($self->id() !~ m/MGRast/) {
1207                            $self->update_stats_for_gap_filling($ElapsedTime);
1208                    }
1209            } else {
1210                    $self->update_model_stats();
1211            }
1212            #Printing the updated SBML file
1213            $self->PrintSBMLFile();
1214            $self->PrintModelLPFile($self->id());
1215            $self->set_status(1,"Auto completion successfully finished");
1216            $self->run_default_model_predictions();
1217    
1218            return $self->success();
1219    }
1220    
1221    =head3 datagapfill
1222    Definition:
1223            success()/fail() = FIGMODELmodel->datagapfill();
1224    Description:
1225            Run gapfilling on the input run specifications
1226    =cut
1227    sub datagapfill {
1228            my ($self,$GapFillingRunSpecs,$TansferFileSuffix) = @_;
1229            my $UniqueFilename = $self->figmodel()->filename();
1230            if (defined($GapFillingRunSpecs) && @{$GapFillingRunSpecs} > 0) {
1231                    print "Gap filling specs:\n".join("\n",@{$GapFillingRunSpecs})."\n\n";
1232                    system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id().$self->selected_version(),"NoBounds",["DataGapFilling"],{"Reactions to knockout" => $self->config("permanently knocked out reactions")->[0],"Gap filling runs" => join(";",@{$GapFillingRunSpecs})},"GapFilling-".$self->id().$self->selected_version()."-".$UniqueFilename.".log",undef,undef));
1233    
1234                    #Checking that the solution exists
1235                    if (!-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt") {
1236                            $self->figmodel()->error_message("FIGMODEL:GapFillingAlgorithm: Could not find MFA output file!");
1237                            $self->figmodel()->database()->print_array_to_file($self->directory().$self->id().$self->selected_version()."-GFS.txt",["Experiment;Solution index;Solution cost;Solution reactions"]);
1238                            return undef;
1239                    }
1240                    my $GapFillResultTable = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt",";","",0,undef);
1241                    if (defined($TansferFileSuffix)) {
1242                            system("cp ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt ".$self->directory().$self->id().$self->selected_version()."-".$TansferFileSuffix.".txt");
1243                            system("cp ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingReport.txt ".$self->directory().$self->id().$self->selected_version()."-".$TansferFileSuffix.".txt");
1244                    }
1245    
1246                    #If the system is not configured to preserve all logfiles, then the mfatoolkit output folder is deleted
1247                    $self->figmodel()->clearing_output($UniqueFilename,"GapFilling-".$self->id().$self->selected_version()."-".$UniqueFilename.".log");
1248                    return $GapFillingRunSpecs;
1249            }
1250            if (defined($TansferFileSuffix)) {
1251                    $self->figmodel()->database()->print_array_to_file($self->directory().$self->id().$self->selected_version()."-".$TansferFileSuffix.".txt",["Experiment;Solution index;Solution cost;Solution reactions"]);
1252            }
1253            return undef;
1254    }
1255    
1256    =head3 datagapgen
1257    Definition:
1258            $model->datagapgen($NumberOfProcessors,$ProcessorIndex,$Filename);
1259    =cut
1260    
1261    sub datagapgen {
1262            my ($self,$Media,$KOList,$NoKOList,$suffix) = @_;
1263            #Setting default values for inputs
1264            if (!defined($NoKOList)) {
1265                    $NoKOList = "none";
1266            }
1267            if (!defined($KOList)) {
1268                    $KOList = "none";
1269            }
1270            #Getting unique filename
1271            my $Filename = $self->figmodel()->filename();
1272            if (!defined($suffix)) {
1273                    $suffix = "-".$Media."-".$KOList."-S";
1274            }
1275            $KOList =~ s/,/;/g;
1276            $NoKOList =~ s/,/;/g;
1277            #Running the gap generation
1278            system($self->figmodel()->GenerateMFAToolkitCommandLineCall($Filename,$self->id().$self->selected_version(),$Media,["GapGeneration"],{"Reactions that should always be active" => $NoKOList,"Reactions to knockout" => $KOList,"Reactions that are always blocked" => "none"},"Gapgeneration-".$self->id().$self->selected_version()."-".$Filename.".log",undef,undef));
1279            my $ProblemReport = $self->figmodel()->LoadProblemReport($Filename);
1280            if (!defined($ProblemReport)) {
1281                    $self->figmodel()->error_message("FIGMODEL:GapGenerationAlgorithm;No problem report;".$Filename.";".$self->id().$self->selected_version().";".$Media.";".$KOList.";".$NoKOList);
1282                    return undef;
1283            }
1284            #Clearing the output folder and log file
1285            $self->figmodel()->clearing_output($Filename,$self->directory()."Gapgeneration-".$self->id().$self->selected_version()."-".$Filename.".log");
1286            #Scheduling the testing of this gap filling solution
1287            my $Tbl = FIGMODELTable->new(["Experiment","Solution index","Solution cost","Solution reactions"],$self->directory().$self->id().$self->selected_version()."-GG".$suffix.".txt",undef,";","|",undef);
1288            for (my $j=0; $j < $ProblemReport->size(); $j++) {
1289                    if ($ProblemReport->get_row($j)->{"Notes"}->[0] =~ m/^Recursive\sMILP\s([^)]+)/) {
1290                            my @SolutionList = split(/\|/,$1);
1291                            for (my $k=0; $k < @SolutionList; $k++) {
1292                                    if ($SolutionList[$k] =~ m/(\d+):(.+)/) {
1293                                            $Tbl->add_row({"Experiment"=>[$Media],"Solution index"=>[$k],"Solution cost"=>[$1],"Solution reactions"=>[$2]});
1294                                    }
1295                            }
1296                    }
1297            }
1298            return $Tbl;
1299    }
1300    
1301    =head3 TestSolutions
1302    Definition:
1303            $model->TestSolutions($ModelID,$NumProcessors,$ProcessorIndex,$GapFill);
1304    Description:
1305    Example:
1306    =cut
1307    
1308    sub TestSolutions {
1309            my ($self,$OriginalErrorFilename,$GapFillResultTable) = @_;
1310            #Getting the filename
1311            my $UniqueFilename = $self->figmodel()->filename();
1312            #Reading in the original error matrix which has the headings for the original model simulation
1313            my $OriginalErrorData;
1314            if (!defined($OriginalErrorFilename) || !-e $self->directory().$OriginalErrorFilename) {
1315                    my ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector) = $self->RunAllStudiesWithDataFast("All");
1316                    $OriginalErrorData = [$HeadingVector,$Errorvector];
1317            } else {
1318                    $OriginalErrorData = $self->figmodel()->database()->load_single_column_file($self->directory().$OriginalErrorFilename,"");
1319            }
1320            my $HeadingHash;
1321            my @HeadingArray = split(/;/,$OriginalErrorData->[0]);
1322            my @OrigErrorArray = split(/;/,$OriginalErrorData->[1]);
1323            for (my $i=0; $i < @HeadingArray; $i++) {
1324                    my @SubArray = split(/:/,$HeadingArray[$i]);
1325                    $HeadingHash->{$SubArray[0].":".$SubArray[1].":".$SubArray[2]} = $i;
1326            }
1327            #Scanning through the gap filling solutions
1328            my $TempVersion = "V".$UniqueFilename;
1329            my $ErrorMatrixLines;
1330            for (my $i=0; $i < $GapFillResultTable->size(); $i++) {
1331                    print "Starting problem solving ".$i."\n";
1332                    my $ErrorLine = $GapFillResultTable->get_row($i)->{"Experiment"}->[0].";".$GapFillResultTable->get_row($i)->{"Solution index"}->[0].";".$GapFillResultTable->get_row($i)->{"Solution cost"}->[0].";".$GapFillResultTable->get_row($i)->{"Solution reactions"}->[0];
1333                    #Integrating solution into test model
1334                    my $ReactionArray;
1335                    my $DirectionArray;
1336                    my @ReactionList = split(/,/,$GapFillResultTable->get_row($i)->{"Solution reactions"}->[0]);
1337                    my %SolutionHash;
1338                    for (my $k=0; $k < @ReactionList; $k++) {
1339                            if ($ReactionList[$k] =~ m/(.+)(rxn\d\d\d\d\d)/) {
1340                                    my $Reaction = $2;
1341                                    my $Sign = $1;
1342                                    if (defined($SolutionHash{$Reaction})) {
1343                                            $SolutionHash{$Reaction} = "<=>";
1344                                    } elsif ($Sign eq "-") {
1345                                            $SolutionHash{$Reaction} = "<=";
1346                                    } elsif ($Sign eq "+") {
1347                                            $SolutionHash{$Reaction} = "=>";
1348                                    } else {
1349                                            $SolutionHash{$Reaction} = $Sign;
1350                                    }
1351                            }
1352                    }
1353                    @ReactionList = keys(%SolutionHash);
1354                    for (my $k=0; $k < @ReactionList; $k++) {
1355                            push(@{$ReactionArray},$ReactionList[$k]);
1356                            push(@{$DirectionArray},$SolutionHash{$ReactionList[$k]});
1357                    }
1358                    print "Integrating solution!\n";
1359                    $self->IntegrateGrowMatchSolution($self->id().$self->selected_version(),$self->directory().$self->id().$TempVersion.".txt",$ReactionArray,$DirectionArray,"Gapfilling ".$GapFillResultTable->get_row($i)->{"Experiment"}->[0],1,1);
1360                    my $model = $self->get_model($self->id().$TempVersion);
1361                    $model->PrintModelLPFile();
1362                    #Running the model against all available experimental data
1363                    print "Running test model!\n";
1364                    my ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector) = $self->RunAllStudiesWithDataFast("All");
1365    
1366                    @HeadingArray = split(/;/,$HeadingVector);
1367                    my @ErrorArray = @OrigErrorArray;
1368                    my @TempArray = split(/;/,$Errorvector);
1369                    for (my $j=0; $j < @HeadingArray; $j++) {
1370                            my @SubArray = split(/:/,$HeadingArray[$j]);
1371                            $ErrorArray[$HeadingHash->{$SubArray[0].":".$SubArray[1].":".$SubArray[2]}] = $TempArray[$j];
1372                    }
1373                    $ErrorLine .= ";".$FalsePostives."/".$FalseNegatives.";".join(";",@ErrorArray);
1374                    push(@{$ErrorMatrixLines},$ErrorLine);
1375                    print "Finishing problem solving ".$i."\n";
1376            }
1377            #Clearing out the test model
1378            if (-e $self->directory().$self->id().$TempVersion.".txt") {
1379                    unlink($self->directory().$self->id().$TempVersion.".txt");
1380                    unlink($self->directory()."SimulationOutput".$self->id().$TempVersion.".txt");
1381            }
1382            return $ErrorMatrixLines;
1383    }
1384    
1385    =head3 status
1386    Definition:
1387            int::model status = FIGMODELmodel->status();
1388    Description:
1389            Returns the current status of the SEED model associated with the input genome ID.
1390            model status = 1: model exists
1391            model status = 0: model is being built
1392            model status = -1: model does not exist
1393            model status = -2: model build failed
1394    =cut
1395    sub status {
1396            my ($self) = @_;
1397            return $self->{_data}->{status}->[0];
1398    }
1399    
1400    =head3 message
1401    Definition:
1402            string::model message = FIGMODELmodel->message();
1403    Description:
1404            Returns a message associated with the models current status
1405    =cut
1406    sub message {
1407            my ($self) = @_;
1408            return $self->{_data}->{message}->[0];
1409    }
1410    
1411    =head3 set_status
1412    Definition:
1413            (success/fail) = FIGMODELmodel->set_status(int::new status,string::status message);
1414    Description:
1415            Changes the current status of the SEED model
1416            new status = 1: model exists
1417            new status = 0: model is being built
1418            new status = -1: model does not exist
1419            new status = -2: model build failed
1420    =cut
1421    sub set_status {
1422            my ($self,$NewStatus,$Message) = @_;
1423    
1424            #Getting the model row from the MODELS table
1425            $self->{_data}->{status}->[0] = $NewStatus;
1426            $self->{_data}->{message}->[0] = $Message;
1427            $self->figmodel()->database()->update_row("MODELS",$self->{_data},"id");
1428            return $self->config("SUCCESS")->[0];
1429    }
1430    
1431    =head3 CreateSingleGenomeReactionList
1432    Definition:
1433            FIGMODEL->CreateSingleGenomeReactionList();
1434    Description:
1435            This function uses fig calls to obtain a list of genes and functions for a genome, and it uses a file mapping reactions and functional roles to produce a reaction list.
1436    Example:
1437    =cut
1438    
1439    sub CreateSingleGenomeReactionList {
1440            my ($self,$RunGapFilling) = @_;
1441    
1442            #Getting genome stats
1443            my $genomestats = $self->figmodel()->get_genome_stats($self->genome());
1444            my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome());
1445            if (!defined($FeatureTable)) {
1446                    $self->figmodel()->error_message("FIGMODEL:CreateSingleGenomeReactionList: ".$self->id()." genome features could not be accessed!");
1447                    return $self->fail();
1448            }
1449            #Checking that the number of genes exceeds the minimum size
1450            if ($FeatureTable->size() < $self->config("minimum genome size for modeling")->[0]) {
1451                    $self->figmodel()->error_message("FIGMODEL:CreateSingleGenomeReactionList: ".$self->id()." genome rejected as too small for modeling!");
1452                    return $self->fail();
1453            }
1454            #Setting up needed variables
1455            my $OriginalModelTable = undef;
1456            #Locking model status table
1457            my $ModelTable = $self->figmodel()->database()->LockDBTable("MODELS");
1458            my $Row = $ModelTable->get_row_by_key($self->id(),"id");
1459            if (!defined($Row) || !defined($Row->{status}->[0])) {
1460                    $self->figmodel()->database()->UnlockDBTable("MODELS");
1461                    $self->figmodel()->error_message("FIGMODEL:CreateSingleGenomeReactionList: ".$self->id()." has no row in models table!");
1462                    return $self->fail();
1463            } elsif ($Row->{status}->[0] == 0) {
1464                    $self->figmodel()->error_message("FIGMODEL:CreateSingleGenomeReactionList:Model is already being built. Canceling current build.");
1465                    $self->figmodel()->database()->UnlockDBTable("MODELS");
1466                    return $self->fail();
1467            }elsif ($Row->{status}->[0] == 1) {
1468                    $OriginalModelTable = $self->reaction_table();
1469                    $self->ArchiveModel();
1470                    $Row->{status}->[0] = 0;
1471                    $Row->{message}->[0] = "Rebuilding preliminary reconstruction";
1472            } else {
1473                    $Row->{status}->[0] = 0;
1474                    $Row->{message}->[0] = "Preliminary reconstruction";
1475            }
1476            #Updating the status table
1477            $self->figmodel()->database()->save_table($ModelTable);
1478            $self->figmodel()->database()->UnlockDBTable("MODELS");
1479            #Sorting GenomeData by gene location on the chromosome
1480            $FeatureTable->sort_rows("MIN LOCATION");
1481            my ($ComplexHash,$SuggestedMappings,$UnrecognizedReactions,$ReactionHash);
1482            my %LastGenePosition;
1483            my $GeneRoles;
1484            for (my $j=0; $j < $FeatureTable->size(); $j++) {
1485                    my $CurrentRow = $FeatureTable->get_row($j);
1486                    #"ID","ALIASES","MIN LOCATION","MAX LOCATION","ROLES","SUBSYSTEMS","SUBSYSTEM CLASS"
1487                    if (defined($CurrentRow)) {
1488                            my $GeneID = $CurrentRow->{"ID"}->[0];
1489                            if ($GeneID =~ m/(peg\.\d+)/) {
1490                                    $GeneID = $1;
1491                            }
1492                            foreach my $Role (@{$CurrentRow->{"ROLES"}}) {
1493                                    if ($self->figmodel()->role_is_valid($Role) != 0) {
1494                                            push(@{$GeneRoles->{$GeneID}},$Role);
1495                                            my $ReactionHashArrayRef = $self->figmodel()->reactions_of_role($Role);
1496                                            if ($ReactionHashArrayRef != 0) {
1497                                                    foreach my $Mapping (@{$ReactionHashArrayRef}) {
1498                                                            if (defined($Mapping->{"REACTION"}) && defined($Mapping->{"MASTER"}) && defined($Mapping->{"SUBSYSTEM"}) && defined($Mapping->{"SOURCE"})) {
1499                                                                    if ($Mapping->{"REACTION"}->[0] =~ m/rxn\d\d\d\d\d/) {
1500                                                                            if ($Mapping->{"MASTER"}->[0] eq 1) {
1501                                                                                    #Creating a complex if consecutive genes have been assigned to the same reaction
1502                                                                                    $ComplexHash->{$Mapping->{"REACTION"}->[0]}->{$Mapping->{"COMPLEX"}->[0]}->{$Role}->{$GeneID} = 1;
1503                                                                                    if (!defined($LastGenePosition{$Mapping->{"REACTION"}->[0]})) {
1504                                                                                            $LastGenePosition{$Mapping->{"REACTION"}->[0]} = $j;
1505                                                                                            push(@{$ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"GENES"}},$GeneID);
1506                                                                                    } elsif (($j-$LastGenePosition{$Mapping->{"REACTION"}->[0]}) < 3 && $LastGenePosition{$Mapping->{"REACTION"}->[0]} != $j) {
1507                                                                                            my $CurrentComplex = pop(@{$ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"GENES"}});
1508                                                                                            push(@{$ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"GENES"}},$CurrentComplex."+".$GeneID);
1509                                                                                    } elsif ($LastGenePosition{$Mapping->{"REACTION"}->[0]} != $j) {
1510                                                                                            push(@{$ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"GENES"}},$GeneID);
1511                                                                                    }
1512                                                                                    $LastGenePosition{$Mapping->{"REACTION"}->[0]} = $j;
1513                                                                                    #Adding a subsystem for the reaction
1514                                                                                    if ($self->figmodel()->subsystem_is_valid($Mapping->{"SUBSYSTEM"}->[0]) == 1) {
1515                                                                                            ($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"SUBSYSTEMS"},my $NumMatches) = $self->figmodel()->add_elements_unique($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"SUBSYSTEMS"},$Mapping->{"SUBSYSTEM"}->[0]);
1516                                                                                            if (!defined($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}) || $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] > 1) {
1517                                                                                                    if ($Mapping->{"SOURCE"}->[0] =~ m/Hope\sFiles/) {
1518                                                                                                            $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] = 1;
1519                                                                                                    } elsif ($Mapping->{"SOURCE"}->[0] =~ m/SEED/) {
1520                                                                                                            $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] = 2;
1521                                                                                                    } elsif (!defined($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}) || $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] > 2) {
1522                                                                                                            $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] = 3;
1523                                                                                                    }
1524                                                                                            }
1525                                                                                    }
1526                                                                                    #Handling confidence
1527                                                                                    if (!defined($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}) || $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] > 2) {
1528                                                                                            if ($Mapping->{"SOURCE"}->[0] =~ m/MATT/) {
1529                                                                                                    $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] = 3;
1530                                                                                            } elsif ($Mapping->{"SOURCE"}->[0] =~ m/CHRIS/) {
1531                                                                                                    $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] = 4;
1532                                                                                            } else {
1533                                                                                                    $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] = 5;
1534                                                                                            }
1535                                                                                    }
1536                                                                                    #Parsing sources
1537                                                                                    ($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"SOURCE"},my $NumMatches) = $self->figmodel()->add_elements_unique($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"SOURCE"},split(/\|/,$Mapping->{"SOURCE"}->[0]));
1538                                                                            } else {
1539                                                                                    push(@{$SuggestedMappings},$GeneID."\t".$Mapping->{"REACTION"}->[0]."\t".$Role);
1540                                                                            }
1541                                                                    } else {
1542                                                                            push(@{$UnrecognizedReactions},$GeneID."\t".$Mapping->{"REACTION"}->[0]."\t".$Role);
1543                                                                    }
1544                                                            }
1545                                                    }
1546                                            }
1547                                    }
1548                            }
1549                    }
1550            }
1551    
1552            #Creating nonadjacent complexes
1553            my @ReactionList = keys(%{$ReactionHash});
1554            foreach my $Reaction (@ReactionList) {
1555                    #If multiple genes are assigned to the reaction, we check if they should should be in a complex
1556                    if (@{$ReactionHash->{$Reaction}->{"GENES"}} > 0 && defined($ComplexHash->{$Reaction})) {
1557                            my $GeneArray;
1558                            foreach my $Complex (keys(%{$ComplexHash->{$Reaction}})) {
1559                                    my %ComplexComponents;
1560                                    foreach my $CurrentGeneSet (@{$ReactionHash->{$Reaction}->{"GENES"}}) {
1561                                            my @GeneList = split(/\+/,$CurrentGeneSet);
1562                                            my %RoleHash;
1563                                            foreach my $Gene (@GeneList) {
1564                                                    foreach my $Role (@{$GeneRoles->{$Gene}}) {
1565                                                            if (defined($ComplexHash->{$Reaction}->{$Complex}->{$Role})) {
1566                                                                    $RoleHash{$Role} = 1;
1567                                                            }
1568                                                    }
1569                                            }
1570                                            if (keys(%RoleHash) > 0) {
1571                                                    if (!defined($ComplexComponents{join("|",sort(keys(%RoleHash)))})) {
1572                                                            my @RoleList = keys(%RoleHash);
1573                                                            my @ComplexList = keys(%ComplexComponents);
1574                                                            foreach my $ComplexSet (@ComplexList) {
1575                                                                    my @RoleList = split(/\|/,$ComplexSet);
1576                                                                    my $Match = 0;
1577                                                                    foreach my $SingleRole (@RoleList) {
1578                                                                            if (defined($RoleHash{$SingleRole})) {
1579                                                                                    $Match = 1;
1580                                                                                    last;
1581                                                                            }
1582                                                                    }
1583                                                                    if ($Match == 1) {
1584                                                                            foreach my $SingleRole (@RoleList) {
1585                                                                                    $RoleHash{$SingleRole} = 1
1586                                                                            }
1587                                                                            push(@{$ComplexComponents{join("|",sort(keys(%RoleHash)))}},@{$ComplexComponents{$ComplexSet}});
1588                                                                            delete $ComplexComponents{$ComplexSet};
1589                                                                    }
1590                                                            }
1591                                                    }
1592                                                    push(@{$ComplexComponents{join("|",sort(keys(%RoleHash)))}},$CurrentGeneSet);
1593                                            }
1594                                    }
1595                                    my @Position;
1596                                    my @Options;
1597                                    my $Count = 0;
1598                                    foreach my $RoleSet (keys(%ComplexComponents)) {
1599                                            push(@Position,0);
1600                                            push(@{$Options[$Count]},@{$ComplexComponents{$RoleSet}});
1601                                            $Count++;
1602                                    }
1603                                    my $Done = 0;
1604                                    $Count = 0;
1605                                    my $NewRelationship;
1606                                    while($Done == 0) {
1607                                            #Creating complex with current indecies
1608                                            $NewRelationship->[$Count] = $Options[0]->[$Position[0]];
1609                                            for (my $i=1; $i < @Position; $i++) {
1610                                                    $NewRelationship->[$Count] .= "+".$Options[$i]->[$Position[$i]];
1611                                            }
1612                                            $NewRelationship->[$Count] = join("+",$self->figmodel()->remove_duplicates(split(/\+/,$NewRelationship->[$Count])));
1613                                            $Count++;
1614                                            #Iterating indecies
1615                                            my $CurrentIndex = 0;
1616                                            while($CurrentIndex >= 0) {
1617                                                    if ($CurrentIndex >= @Position) {
1618                                                            $CurrentIndex = -1000;
1619                                                    } elsif ($Position[$CurrentIndex]+1 == @{$Options[$CurrentIndex]}) {
1620                                                            $Position[$CurrentIndex] = -1;
1621                                                            $CurrentIndex++;
1622                                                    } else {
1623                                                            $Position[$CurrentIndex]++;
1624                                                            $CurrentIndex--;
1625                                                    }
1626                                            }
1627                                            if ($CurrentIndex == -1000) {
1628                                                    $Done = 1;
1629                                            }
1630                                    }
1631                                    push(@{$GeneArray},@{$NewRelationship});
1632                            }
1633                            @{$ReactionHash->{$Reaction}->{"GENES"}} = $self->figmodel()->remove_duplicates(@{$GeneArray});
1634                    }
1635            }
1636    
1637            #Getting the reaction table
1638            my $ReactionTable = $self->figmodel()->database()->GetDBTable("REACTIONS");
1639    
1640            #Creating the model reaction table
1641            my $NewModelTable = FIGMODELTable->new(["LOAD","DIRECTIONALITY","COMPARTMENT","ASSOCIATED PEG","SUBSYSTEM","CONFIDENCE","REFERENCE","NOTES"],$self->directory().$self->id().".txt",["LOAD"],";","|","REACTIONS\n");
1642            @ReactionList = keys(%{$ReactionHash});
1643            foreach my $ReactionID (@ReactionList) {
1644                    #Getting the thermodynamic reversibility from the database
1645                    my $Directionality = $self->figmodel()->reversibility_of_reaction($ReactionID);
1646                    my $Subsystem = "NONE";
1647                    if (defined($ReactionHash->{$ReactionID}->{"SUBSYSTEMS"})) {
1648                            $Subsystem = join("|",@{$ReactionHash->{$ReactionID}->{"SUBSYSTEMS"}});
1649                    }
1650                    my $Source = "NONE";
1651                    if (defined($ReactionHash->{$ReactionID}->{"SOURCE"})) {
1652                            $Source = join("|",@{$ReactionHash->{$ReactionID}->{"SOURCE"}});
1653                    }
1654                    $NewModelTable->add_row({"LOAD" => [$ReactionID],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => [join("|",@{$ReactionHash->{$ReactionID}->{"GENES"}})],"SUBSYSTEM" => [$Subsystem],"CONFIDENCE" => [$ReactionHash->{$ReactionID}->{"CONFIDENCE"}->[0]],"REFERENCE" => [$Source],"NOTES" => ["NONE"]});
1655            }
1656    
1657            #Adding the spontaneous and universal reactions
1658            foreach my $ReactionID (@{$self->config("spontaneous reactions")}) {
1659                    #Getting the thermodynamic reversibility from the database
1660                    my $Directionality = $self->figmodel()->reversibility_of_reaction($ReactionID);
1661                    #Checking if the reaction is already in the model
1662                    if (!defined($NewModelTable->get_row_by_key($ReactionID,"LOAD"))) {
1663                            $NewModelTable->add_row({"LOAD" => [$ReactionID],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["SPONTANEOUS"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [4],"REFERENCE" => ["SPONTANEOUS"],"NOTES" => ["NONE"]});
1664                    }
1665            }
1666            foreach my $ReactionID (@{$self->config("universal reactions")}) {
1667                    #Getting the thermodynamic reversibility from the database
1668                    my $Directionality = $self->figmodel()->reversibility_of_reaction($ReactionID);
1669                    #Checking if the reaction is already in the model
1670                    if (!defined($NewModelTable->get_row_by_key($ReactionID,"LOAD"))) {
1671                            $NewModelTable->add_row({"LOAD" => [$ReactionID],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["UNIVERSAL"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [4],"REFERENCE" => ["UNIVERSAL"],"NOTES" => ["NONE"]});
1672                    }
1673            }
1674    
1675            #Checking if a biomass reaction already exists
1676            my $BiomassReactionRow = $self->figmodel()->database()->get_row_by_key("BIOMASS TABLE",$self->id(),"MODELS");
1677            if (!defined($BiomassReactionRow)) {
1678                    $BiomassReactionRow = $self->BuildSpecificBiomassReaction();
1679                    if (!defined($BiomassReactionRow)) {
1680                            $self->set_status(-2,"Preliminary reconstruction failed: could not generate biomass reaction");
1681                            $self->figmodel()->error_message("FIGMODEL:CreateModelReactionList: Could not generate biomass function for ".$self->id()."!");
1682                            return $self->fail();
1683                    }
1684            }
1685            my $ReactionList = $BiomassReactionRow->{"ESSENTIAL REACTIONS"};
1686            push(@{$ReactionList},$BiomassReactionRow->{DATABASE}->[0]);
1687    
1688            #Adding biomass reactions to the model table
1689            foreach my $BOFReaction (@{$ReactionList}) {
1690                    #Getting the thermodynamic reversibility from the database
1691                    my $Directionality = $self->figmodel()->reversibility_of_reaction($BOFReaction);
1692                    #Checking if the reaction is already in the model
1693                    if (!defined($NewModelTable->get_row_by_key($BOFReaction,"LOAD"))) {
1694                            if ($BOFReaction =~ m/bio/) {
1695                                    $NewModelTable->add_row({"LOAD" => [$BOFReaction],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["BOF"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [1],"REFERENCE" => ["Biomass objective function"],"NOTES" => ["NONE"]});
1696                            } else {
1697                                    $NewModelTable->add_row({"LOAD" => [$BOFReaction],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["INITIAL GAP FILLING"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [5],"REFERENCE" => ["Initial gap filling"],"NOTES" => ["NONE"]});
1698                            }
1699                    }
1700            }
1701    
1702            #Completing any incomplete reactions sets
1703            my $ReactionSetTable = $self->figmodel()->database()->GetDBTable("REACTION SETS");
1704            for (my $i=0; $i < $ReactionSetTable->size(); $i++) {
1705                    if (defined($NewModelTable->get_row_by_key($ReactionSetTable->get_row($i)->{"Trigger reaction"}->[0],"LOAD"))) {
1706                            foreach my $Reaction (@{$ReactionSetTable->get_row($i)->{"Dependant reactions"}}) {
1707                                    if (!defined($NewModelTable->get_row_by_key($ReactionSetTable->get_row($i)->{"Trigger reaction"}->[0],"LOAD"))) {
1708                                            #Getting the thermodynamic reversibility from the database
1709                                            my $Directionality = $self->figmodel()->reversibility_of_reaction($Reaction);
1710                                            $NewModelTable->add_row({"LOAD" => [$Reaction],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["REACTION SET GAP FILLING"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [5],"REFERENCE" => ["Added due to presence of ".$ReactionSetTable->get_row($i)->{"Trigger reaction"}->[0]],"NOTES" => ["NONE"]});
1711                                    }
1712                            }
1713                    }
1714            }
1715    
1716            #Now we compare the model to the previous model to determine if any differences exist that aren't gap filling reactions
1717            if (defined($OriginalModelTable)) {
1718                    my $PerfectMatch = 1;
1719                    my $ReactionCount = 0;
1720                    for (my $i=0; $i < $OriginalModelTable->size(); $i++) {
1721                            #We only check that nongapfilling reactions exist in the new model
1722                            if ($OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0] !~ m/GAP/ || $OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0] eq "INITIAL GAP FILLING") {
1723                                    $ReactionCount++;
1724                                    my $Row = $NewModelTable->get_row_by_key($OriginalModelTable->get_row($i)->{"LOAD"}->[0],"LOAD");
1725                                    if (defined($Row)) {
1726                                            #We check that the reaction directionality is identical
1727                                            if ($Row->{"DIRECTIONALITY"}->[0] ne $OriginalModelTable->get_row($i)->{"DIRECTIONALITY"}->[0]) {
1728                                                    if (defined($OriginalModelTable->get_row($i)->{"NOTES"}->[0]) && $OriginalModelTable->get_row($i)->{"NOTES"}->[0] =~ m/Directionality\sswitched\sfrom\s([^\s])/) {
1729                                                            if ($1 ne $Row->{"DIRECTIONALITY"}->[0]) {
1730                                                                    print "Directionality mismatch for reaction ".$OriginalModelTable->get_row($i)->{"LOAD"}->[0].": ".$1." vs ".$Row->{"DIRECTIONALITY"}->[0]."\n";
1731                                                                    $PerfectMatch = 0;
1732                                                                    last;
1733                                                            }
1734                                                    } else {
1735                                                            print "Directionality mismatch for reaction ".$OriginalModelTable->get_row($i)->{"LOAD"}->[0].": ".$OriginalModelTable->get_row($i)->{"DIRECTIONALITY"}->[0]." vs ".$Row->{"DIRECTIONALITY"}->[0]."\n";
1736                                                            $PerfectMatch = 0;
1737                                                            last;
1738                                                    }
1739                                            }
1740                                            #We check that the genes assigned to the reaction are identical
1741                                            if ($PerfectMatch == 1 && @{$OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}} != @{$Row->{"ASSOCIATED PEG"}}) {
1742                                                    print "Gene associatation mismatch for reaction ".$OriginalModelTable->get_row($i)->{"LOAD"}->[0].": ".@{$OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}}." vs ".@{$Row->{"ASSOCIATED PEG"}}."\n";
1743                                                    $PerfectMatch = 0;
1744                                                    last;
1745                                            }
1746                                            if ($PerfectMatch == 1) {
1747                                                    my @GeneSetOne = sort(@{$OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}});
1748                                                    my @GeneSetTwo = sort(@{$Row->{"ASSOCIATED PEG"}});
1749                                                    for (my $j=0; $j < @GeneSetOne; $j++) {
1750                                                            if ($GeneSetOne[$j] ne $GeneSetTwo[$j]) {
1751                                                                    print "Gene mismatch for reaction ".$OriginalModelTable->get_row($i)->{"LOAD"}->[0].": ".$GeneSetOne[$j]." vs ".$GeneSetTwo[$j]."\n";
1752                                                                    $PerfectMatch = 0;
1753                                                                    $i = $OriginalModelTable->size();
1754                                                                    last;
1755                                                            }
1756                                                    }
1757                                            }
1758                                    } else {
1759                                            print "Original model contains an extra reaction:".$OriginalModelTable->get_row($i)->{"LOAD"}->[0]."\n";
1760                                            $PerfectMatch = 0;
1761                                            last;
1762                                    }
1763                            }
1764                    }
1765                    if ($PerfectMatch == 1 && $ReactionCount == $NewModelTable->size()) {
1766                            #Bailing out of function as the model has not changed
1767                            $self->set_status(1,"rebuild canceled because model has not changed");
1768                            return $self->success();
1769                    }
1770            }
1771    
1772            #Saving the new model to file
1773            $self->set_status(1,"Preliminary reconstruction complete");
1774            $self->figmodel()->database()->save_table($NewModelTable);
1775            $self->{_reaction_data} = $NewModelTable;
1776            #Clearing the previous model from the cache
1777            $self->figmodel()->ClearDBModel($self->id(),1);
1778            #Updating the model stats table
1779            $self->update_stats_for_build();
1780            $self->PrintSBMLFile();
1781    
1782            #Adding model to gapfilling queue
1783            if (defined($RunGapFilling) && $RunGapFilling == 1) {
1784                    $self->set_status(1,"Autocompletion queued");
1785                    $self->figmodel()->add_job_to_queue("gapfillmodel?".$self->id(),"QSUB","cplex",$self->owner(),"BACK");
1786            }
1787            return $self->success();
1788    }
1789    
1790    =head3 CreateMetaGenomeReactionList
1791    Definition:
1792            (success/fail) = FIGMODELmodel->CreateMetaGenomeReactionList();
1793    Description:
1794            This is the code called to create or update the reaction list for a metgenome model
1795    =cut
1796    
1797    sub CreateMetaGenomeReactionList {
1798            my ($self) = @_;
1799    
1800            #Checking if the metagenome file exists
1801            if (!-e $self->config("raw MGRAST directory")->[0].$self->genome().".summary") {
1802                    $self->error_message("FIGMODEL:CreateMetaGenomeReactionList: could not find raw data file for metagenome ".$self->genome());
1803            }
1804            #Loading metagenome data
1805            my $MGRASTData = $self->figmodel()->database()->load_multiple_column_file($self->config("raw MGRAST directory")->[0].$self->genome().".summary","\t");
1806            if (!defined($MGRASTData)) {
1807                    $self->error_message("FIGMODEL:CreateMetaGenomeReactionList: could not find raw data file for metagenome ".$self->genome());
1808            }
1809    
1810            #Setting up needed variables
1811            my $OriginalModelTable = undef;
1812    
1813            #Checking status
1814            if ($self->status() < 0) {
1815                    $self->set_status(0,"Preliminary reconstruction");
1816            } elsif ($self->status() == 0) {
1817                    $self->error_message("FIGMODEL->CreateModelReactionList:Model is already being built. Canceling current build.");
1818                    return $self->fail();
1819            } else {
1820                    $OriginalModelTable = $self->reaction_table();
1821                    $self->ArchiveModel();
1822                    $self->set_status(0,"Rebuilding preliminary reconstruction");
1823            }
1824    
1825            #Getting the reaction table
1826            my $ReactionTable = $self->figmodel()->database()->GetDBTable("REACTIONS");
1827            #Creating model table
1828            my $ModelTable = FIGMODELTable->new(["LOAD","DIRECTIONALITY","COMPARTMENT","ASSOCIATED PEG","SUBSYSTEM","CONFIDENCE","REFERENCE","NOTES"],$self->directory().$self->id().".txt",["LOAD"],";","|","REACTIONS\n");
1829            for (my $i=0; $i < @{$MGRASTData};$i++) {
1830                    #MD5,PEG,number of sims,role,sim e-scores
1831                    my $Role = $MGRASTData->[$i]->[3];
1832                    my $MD5 = $MGRASTData->[$i]->[0];
1833                    my $peg = $MGRASTData->[$i]->[1];
1834                    my $sims = $MGRASTData->[$i]->[4];
1835                    $sims =~ s/;/,/g;
1836    
1837                    #Checking for subsystems
1838                    my $GeneSubsystems = $self->figmodel()->subsystems_of_role($Role);
1839                    #Checking if there are reactions associated with this role
1840                    my $ReactionHashArrayRef = $self->figmodel()->reactions_of_role($Role);
1841                    if ($ReactionHashArrayRef != 0) {
1842                            foreach my $Mapping (@{$ReactionHashArrayRef}) {
1843                                    if (defined($Mapping->{"REACTION"}) && defined($Mapping->{"MASTER"}) && defined($Mapping->{"SUBSYSTEM"}) && defined($Mapping->{"SOURCE"})) {
1844                                            if ($Mapping->{"REACTION"}->[0] =~ m/rxn\d\d\d\d\d/) {
1845                                                    if ($Mapping->{"MASTER"}->[0] eq 1) {
1846                                                            #Checking if the reaction is already in the model
1847                                                            my $ReactionRow = $ModelTable->get_row_by_key($Mapping->{"REACTION"}->[0],"LOAD");
1848                                                            if (!defined($ReactionRow)) {
1849                                                                    $ReactionRow = {"LOAD" => [$Mapping->{"REACTION"}->[0]],"DIRECTIONALITY" => [$self->figmodel()->reversibility_of_reaction($Mapping->{"REACTION"}->[0])],"COMPARTMENT" => ["c"]};
1850                                                                    $ModelTable->add_row($ReactionRow);
1851                                                            }
1852                                                            push(@{$ReactionRow->{"ASSOCIATED PEG"}},substr($peg,4));
1853                                                            push(@{$ReactionRow->{"REFERENCE"}},$MD5.":".$Role);
1854                                                            push(@{$ReactionRow->{"CONFIDENCE"}},$sims);
1855                                                            if (defined($GeneSubsystems)) {
1856                                                                    push(@{$ReactionRow->{"SUBSYSTEM"}},@{$GeneSubsystems});
1857                                                            }
1858                                                    }
1859                                            }
1860                                    }
1861                            }
1862                    }
1863            }
1864    
1865            #Adding the spontaneous and universal reactions
1866            foreach my $ReactionID (@{$self->config("spontaneous reactions")}) {
1867                    #Getting the thermodynamic reversibility from the database
1868                    my $Directionality = $self->figmodel()->reversibility_of_reaction($ReactionID);
1869                    #Checking if the reaction is already in the model
1870                    if (!defined($ModelTable->get_row_by_key($ReactionID,"LOAD"))) {
1871                            $ModelTable->add_row({"LOAD" => [$ReactionID],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["SPONTANEOUS"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [4],"REFERENCE" => ["SPONTANEOUS"],"NOTES" => ["NONE"]});
1872                    }
1873            }
1874            foreach my $ReactionID (@{$self->config("universal reactions")}) {
1875                    #Getting the thermodynamic reversibility from the database
1876                    my $Directionality = $self->figmodel()->reversibility_of_reaction($ReactionID);
1877                    #Checking if the reaction is already in the model
1878                    if (!defined($ModelTable->get_row_by_key($ReactionID,"LOAD"))) {
1879                            $ModelTable->add_row({"LOAD" => [$ReactionID],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["UNIVERSAL"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [4],"REFERENCE" => ["UNIVERSAL"],"NOTES" => ["NONE"]});
1880                    }
1881            }
1882    
1883            #Completing any incomplete reactions sets
1884            my $ReactionSetTable = $self->figmodel()->database()->GetDBTable("REACTION SETS");
1885            for (my $i=0; $i < $ReactionSetTable->size(); $i++) {
1886                    if (defined($ModelTable->get_row_by_key($ReactionSetTable->get_row($i)->{"Trigger reaction"}->[0],"LOAD"))) {
1887                            foreach my $Reaction (@{$ReactionSetTable->get_row($i)->{"Dependant reactions"}}) {
1888                                    if (!defined($ModelTable->get_row_by_key($ReactionSetTable->get_row($i)->{"Trigger reaction"}->[0],"LOAD"))) {
1889                                            #Getting the thermodynamic reversibility from the database
1890                                            my $Directionality = $self->figmodel()->reversibility_of_reaction($Reaction);
1891                                            $ModelTable->add_row({"LOAD" => [$Reaction],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["REACTION SET GAP FILLING"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [5],"REFERENCE" => ["Added due to presence of ".$ReactionSetTable->get_row($i)->{"Trigger reaction"}->[0]],"NOTES" => ["NONE"]});
1892                                    }
1893                            }
1894                    }
1895            }
1896    
1897            #Clearing the previous model from the cache
1898            $self->figmodel()->ClearDBModel($self->id(),1);
1899            $ModelTable->save();
1900    
1901            return $self->success();
1902    }
1903    
1904    =head3 ArchiveModel
1905    Definition:
1906            (success/fail) = FIGMODELmodel->ArchiveModel();
1907    Description:
1908            This function archives the specified model in the model directory with the current version numbers appended.
1909            This function is used to preserve old versions of models prior to overwriting so new versions may be compared with old versions.
1910    =cut
1911    sub ArchiveModel {
1912            my ($self) = @_;
1913    
1914            #Making sure the model exists
1915            if (!defined($self->stats())) {
1916                    $self->figmodel()->error_message("FIGMODEL:ArchiveModel: Could not find specified ".$self->id()." in database!");
1917                    return $self->fail();
1918            }
1919    
1920            #Checking that the model file exists
1921            if (!(-e $self->filename())) {
1922                    $self->figmodel()->error_message("FIGMODEL:ArchiveModel: Model file ".$self->filename()." not found!");
1923                    return $self->fail();
1924            }
1925    
1926            #Copying the model file
1927            system("cp ".$self->filename()." ".$self->directory().$self->id().$self->version().".txt");
1928    }
1929    
1930    =head3 PrintModelDataToFile
1931    Definition:
1932            (success/fail) = FIGMODELmodel->PrintModelDataToFile();
1933    Description:
1934            This function uses the MFAToolkit to print out all of the compound and reaction data for the input model.
1935            Some of the data printed by the toolkit is calculated internally in the toolkit and not stored in any files, so this data can only be retrieved through this
1936            function. The LoadModel function for example would not give you this data.
1937    =cut
1938    sub PrintModelDataToFile {
1939            my($self) = @_;
1940    
1941            #Running the MFAToolkit on the model file
1942            my $OutputIndex = $self->figmodel()->filename();
1943            my $Command = $self->config("MFAToolkit executable")->[0]." parameterfile ../Parameters/Printing.txt resetparameter output_folder ".$OutputIndex.'/ LoadCentralSystem "'.$self->filename().'"';
1944            system($Command);
1945    
1946            #Copying the model file printed by the toolkit out of the output directory and into the model directory
1947            if (!-e $self->config("MFAToolkit output directory")->[0].$OutputIndex."/".$self->id().$self->selected_version().".txt") {
1948                    $self->figmodel()->error_message("New model file not created due to an error. Check that the input modelfile exists.");
1949                    $self->figmodel()->cleardirectory($OutputIndex);
1950                    return $self->fail();
1951            }
1952    
1953            $Command = 'cp "'.$self->config("MFAToolkit output directory")->[0].$OutputIndex."/".$self->id().$self->selected_version().'.txt" "'.$self->directory().$self->id().$self->selected_version().'Data.txt"';
1954            system($Command);
1955            $Command = 'cp "'.$self->config("MFAToolkit output directory")->[0].$OutputIndex.'/ErrorLog0.txt" "'.$self->directory().'ModelErrors.txt"';
1956            system($Command);
1957            $self->figmodel()->cleardirectory($OutputIndex);
1958            return $self->success();
1959    }
1960    
1961    =head2 Analysis Functions
1962    
1963    =head3 run_microarray_analysis
1964    Definition:
1965            int::status = FIGMODEL->run_microarray_analysis(string::media,string::job id,string::gene calls);
1966    Description:
1967            Runs microarray analysis attempting to turn off genes that are inactive in the microarray
1968    =cut
1969    sub run_microarray_analysis {
1970            my ($self,$media,$jobid,$index,$genecall) = @_;
1971            $genecall =~ s/_/:/g;
1972            $genecall =~ s/\//;/g;
1973            #print "\n\n".$genecall."\n\n";
1974            my $command = $self->figmodel()->GenerateMFAToolkitCommandLineCall($jobid,$self->id(),$media,["MFA","MicroarrayAssertions"],{"Microarray assertions" => $self->id().";".$index.";".$genecall,"MFASolver" => "CPLEX","Network output location" => "/scratch/"},"MicroarrayAnalysis-".$jobid.".txt",undef,$self->selected_version());
1975            #print $command."\n";
1976            system($command);
1977            #system("/home/chenry/Software/MFAToolkitRepository/Linux/mfatoolkit resetparameter \"user bounds filename\" \"Carbon-D-Glucose.txt\" resetparameter output_folder \"32749.354149.0/\" resetparameter \"Microarray assertions\" \"Seed83333.1;0;peg.3130:-1;peg.4035:-1\" resetparameter \"MFASolver\" \"CPLEX\" resetparameter \"Network output location\" \"/scratch/\" LoadCentralSystem \"/vol/model-dev/MODEL_DEV_DB/Models/83333.1/Seed83333.1.txt\" > \"/vol/model-dev/MODEL_DEV_DB/ReactionDB/log/MicroarrayAnalysis-32749.354149.0.txt\"");
1978    }
1979    
1980    =head3 find_minimal_pathways
1981    Definition:
1982            int::status = FIGMODEL->find_minimal_pathways(string::media,string::objective);
1983    Description:
1984            Runs microarray analysis attempting to turn off genes that are inactive in the microarray
1985    =cut
1986    sub find_minimal_pathways {
1987            my ($self,$media,$objective,$solutionnum) = @_;
1988    
1989            #Setting default media
1990            if (!defined($media)) {
1991                    $media = "Complete";
1992            }
1993    
1994            #Setting default solution number
1995            if (!defined($solutionnum)) {
1996                    $solutionnum = "5";
1997            }
1998    
1999            #Translating objective
2000            my $objectivestring;
2001            if ($objective eq "ALL") {
2002                    #Getting the list of universal building blocks
2003                    my $buildingblocks = $self->config("universal building blocks");
2004                    my @objectives = keys(%{$buildingblocks});
2005                    #Getting the nonuniversal building blocks
2006                    my $otherbuildingblocks = $self->config("nonuniversal building blocks");
2007                    my @array = keys(%{$otherbuildingblocks});
2008                    if (defined($self->get_biomass()) && defined($self->figmodel()->get_reaction($self->get_biomass()->{"LOAD"}->[0]))) {
2009                            my $equation = $self->figmodel()->get_reaction($self->get_biomass()->{"LOAD"}->[0])->{"EQUATION"}->[0];
2010                            if (defined($equation)) {
2011                                    for (my $i=0; $i < @array; $i++) {
2012                                            if (CORE::index($equation,$array[$i]) > 0) {
2013                                                    push(@objectives,$array[$i]);
2014                                            }
2015                                    }
2016                            }
2017                    }
2018                    for (my $i=0; $i < @objectives; $i++) {
2019                            $self->find_minimal_pathways($media,$objectives[$i])
2020                    }
2021            } elsif ($objective eq "ENERGY") {
2022                    $objectivestring = "MAX;FLUX;rxn00062;c;1";
2023            } elsif ($objective =~ m/cpd\d\d\d\d\d/) {
2024                    $objectivestring = "MIN;DRAIN_FLUX;".$objective.";c;1";
2025            }
2026    
2027            #Setting additional exchange fluxes
2028            my $additionalexchange = $self->config("default exchange fluxes")->[0];
2029            if ($objective eq "cpd15665" || $objective eq "cpd15667" || $objective eq "cpd15668" || $objective eq "cpd15669") {
2030                    $additionalexchange .= ";cpd15666[c]:0:100";
2031            } elsif ($objective eq "cpd11493") {
2032                    $additionalexchange .= ";cpd12370[c]:0:100";
2033            } elsif ($objective eq "cpd00166") {
2034                    $additionalexchange .= ";cpd01997[c]:0:100;cpd03422[c]:0:100";
2035            }
2036    
2037            #Running MFAToolkit
2038            my $filename = $self->figmodel()->filename();
2039            my $command = $self->figmodel()->GenerateMFAToolkitCommandLineCall($filename,$self->id(),$media,["MFA"],{"Recursive MILP solution limit" => $solutionnum,"Generate pathways to objective" => 1,"MFASolver" => "CPLEX","objective" => $objectivestring,"exchange species" => $additionalexchange},"MinimalPathways-".$media."-".$self->id().$self->selected_version().".txt",undef,$self->selected_version());
2040            system($command);
2041    
2042            #Loading problem report
2043            my $results = $self->figmodel()->LoadProblemReport($filename);
2044            #Clearing output
2045            $self->figmodel()->clearing_output($filename,"MinimalPathways-".$media."-".$self->id()."-".$objective.".txt");
2046            if (!defined($results)) {
2047                    return;
2048            }
2049    
2050            #Parsing output
2051            my @Array;
2052            my $row = $results->get_row(1);
2053            if (defined($row->{"Notes"}->[0])) {
2054                    $_ = $row->{"Notes"}->[0];
2055                    @Array = /\d+:([^\|]+)\|/g;
2056            }
2057    
2058            # Storing data in a figmodel table
2059            my $TableObject;
2060            if (-e $self->directory()."MinimalPathways-".$media."-".$self->id().$self->selected_version().".txt") {
2061                    $TableObject->load_table($self->directory()."MinimalPathways-".$media."-".$self->id().$self->selected_version().".txt",";","|",0,["OBJECTIVE"]);
2062            } else {
2063                    $TableObject = FIGMODELTable->new(["OBJECTIVE","REACTIONS"],$self->directory()."MinimalPathways-".$media."-".$self->id().$self->selected_version().".txt",["OBJECTIVE"],";","|",undef);
2064            }
2065            my $tablerow = $TableObject->get_row_by_key($objective,"OBJECTIVE",1);
2066            push(@{$tablerow->{"REACTIONS"}},@Array);
2067            $TableObject->save();
2068    }
2069    
2070    =head3 calculate_growth
2071    Definition:
2072            string::growth = FIGMODELmodel->calculate_growth(string:media);
2073    Description:
2074            Calculating growth in the input media
2075    =cut
2076    sub calculate_growth {
2077            my ($self,$Media) = @_;
2078            my $UniqueFilename = $self->figmodel()->filename();
2079            system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$Media,["ProductionMFA"],{"optimize metabolite production if objective is zero" => 1},$self->id()."-".$Media."-GrowthTest.txt",undef,$self->selected_version()));
2080            my $ProblemReport = $self->figmodel()->LoadProblemReport($UniqueFilename);
2081            my $Result;
2082            if (defined($ProblemReport)) {
2083                    my $Row = $ProblemReport->get_row(0);
2084                    if (defined($Row) && defined($Row->{"Objective"}->[0])) {
2085                            if ($Row->{"Objective"}->[0] < 0.00000001) {
2086                                    $Result = "NOGROWTH:".$Row->{"Individual metabolites with zero production"}->[0];
2087                            } else {
2088                                    $Result = $Row->{"Objective"}->[0];
2089                            }
2090                    }
2091            }
2092            return $Result;
2093    }
2094    
2095    =head3 classify_model_reactions
2096    Definition:
2097            (FIGMODELTable:Reaction classes,FIGMODELTable:Compound classes) = FIGMODELmodel->classify_model_reactions(string:media);
2098    Description:
2099            This function uses the MFAToolkit to minimize and maximize the flux through every reaction in the input model during minimal growth on the input media.
2100            The results are returned in a hash of strings where the keys are the reaction IDs and the strings are structured as follows: "Class;Min flux;Max flux".
2101            Possible values for "Class" include:
2102            1.) Positive: these reactions are essential in the forward direction.
2103            2.) Negative: these reactions are essential in the reverse direction.
2104            3.) Positive variable: these reactions are nonessential, but they only ever proceed in the forward direction.
2105            4.) Negative variable: these reactions are nonessential, but they only ever proceed in the reverse direction.
2106            5.) Variable: these reactions are nonessential and proceed in the forward or reverse direction.
2107            6.) Blocked: these reactions never carry any flux at all in the media condition tested.
2108            7.) Dead: these reactions are disconnected from the network.
2109    =cut
2110    sub classify_model_reactions {
2111            my ($self,$Media,$SaveChanges) = @_;
2112    
2113            #Getting unique file for printing model output
2114            my $UniqueFilename = $self->figmodel()->filename();
2115            #Running the MFAToolkit
2116            system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$Media,["ProductionMFA"],{"identify dead ends" => 1,"find tight bounds" => 1,"MFASolver" => "GLPK"},"Classify-".$self->id().$self->selected_version()."-".$UniqueFilename.".log",undef,$self->selected_version()));
2117            #Reading in the output bounds file
2118            my ($ReactionTB,$CompoundTB,$DeadCompounds,$DeadEndCompounds,$DeadReactions);
2119            if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsReactionData0.txt") {
2120                    $ReactionTB = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsReactionData0.txt",";","|",1,["DATABASE ID"]);
2121            }
2122            if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsCompoundData0.txt") {
2123                    $CompoundTB = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsCompoundData0.txt",";","|",1,["DATABASE ID"]);
2124            }
2125            if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/DeadReactions.txt") {
2126                    $DeadReactions = $self->figmodel()->put_array_in_hash($self->figmodel()->database()->load_single_column_file($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/DeadReactions.txt",""));
2127            }
2128            if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/DeadMetabolites.txt") {
2129                    $DeadCompounds = $self->figmodel()->put_array_in_hash($self->figmodel()->database()->load_single_column_file($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/DeadMetabolites.txt",""));
2130            }
2131            if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/DeadEndMetabolites.txt") {
2132                    $DeadEndCompounds = $self->figmodel()->put_array_in_hash($self->figmodel()->database()->load_single_column_file($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/DeadEndMetabolites.txt",""));
2133            }
2134            if (!defined($ReactionTB) && !defined($CompoundTB)) {
2135                    print STDERR "FIGMODEL:ClassifyModelReactions: Classification file not found when classifying reactions in ".$self->id().$self->selected_version()." with ".$Media." media. Most likely the model did not grow.\n";
2136                    return (undef,undef);
2137            }
2138    
2139            #Clearing output
2140            $self->figmodel()->clearing_output($UniqueFilename,"Classify-".$self->id().$self->selected_version()."-".$UniqueFilename.".log");
2141            #Creating the table objects that will hold the results of the reaction classification
2142            my $rxnclasstable = $self->reaction_class_table();
2143            my $cpdclasstable = $self->compound_class_table();
2144            #Loading the compound table
2145            if (defined($CompoundTB)) {
2146                    for (my $i=0; $i < $CompoundTB->size(); $i++) {
2147                            my $Row = $CompoundTB->get_row($i);
2148                            if (defined($Row->{"DATABASE ID"})) {
2149                                    #Getting the compound row
2150                                    my $CpdRow = $cpdclasstable->get_row_by_key($Row->{"DATABASE ID"}->[0].$Row->{COMPARTMENT}->[0],"COMPOUND",1);
2151                                    #Setting row values
2152                                    my $Max = 0;
2153                                    my $Min = 0;
2154                                    my $Class = "Unknown";
2155                                    if (defined($DeadCompounds) && defined($DeadCompounds->{$Row->{"DATABASE ID"}->[0]})) {
2156                                            $Class = "Dead";
2157                                    } elsif (defined($DeadEndCompounds) && defined($DeadEndCompounds->{$Row->{"DATABASE ID"}->[0]})) {
2158                                            $Class = "Deadend";
2159                                    } elsif (defined($Row->{"Min DRAIN_FLUX"}) && defined($Row->{"Max DRAIN_FLUX"}) && $Row->{"Min DRAIN_FLUX"}->[0] ne "1e+07") {
2160                                            $Max = $Row->{"Max DRAIN_FLUX"}->[0];
2161                                            $Min = $Row->{"Min DRAIN_FLUX"}->[0];
2162                                            if ($Row->{"Min DRAIN_FLUX"}->[0] > 0.00000001) {
2163                                                    $Class = "Positive";
2164                                            } elsif ($Row->{"Max DRAIN_FLUX"}->[0] < -0.00000001) {
2165                                                    $Class = "Negative";
2166                                            } elsif ($Row->{"Min DRAIN_FLUX"}->[0] < -0.00000001) {
2167                                                    if ($Row->{"Max DRAIN_FLUX"}->[0] > 0.00000001) {
2168                                                            $Class = "Variable";
2169                                                    } else {
2170                                                            $Max = 0;
2171                                                            $Class = "Negative variable";
2172                                                    }
2173                                            } elsif ($Row->{"Max DRAIN_FLUX"}->[0] > 0.00000001) {
2174                                                    $Min = 0;
2175                                                    $Class = "Positive variable";
2176                                            } else {
2177                                                    $Min = 0;
2178                                                    $Max = 0;
2179                                                    $Class = "Blocked";
2180                                            }
2181                                    }
2182                                    my $index = 0;
2183                                    if (defined($CpdRow->{MEDIA})) {
2184                                            for (my $i=0; $i < @{$CpdRow->{MEDIA}};$i++) {
2185                                                    $index++;
2186                                                    if ($CpdRow->{MEDIA}->[$i] eq $Media) {
2187                                                            $index = $i;
2188                                                            last;
2189                                                    }
2190                                            }
2191                                    }
2192                                    $CpdRow->{MIN}->[$index] = $Min;
2193                                    $CpdRow->{MAX}->[$index] = $Max;
2194                                    $CpdRow->{CLASS}->[$index] = $Class;
2195                                    $CpdRow->{MEDIA}->[$index] = $Media;
2196                            }
2197                    }
2198                    if (!defined($SaveChanges) || $SaveChanges == 1) {
2199                            $cpdclasstable->save();
2200                    }
2201            }
2202            if (defined($ReactionTB)) {
2203                    for (my $i=0; $i < $ReactionTB->size(); $i++) {
2204                            my $Row = $ReactionTB->get_row($i);
2205                            if (defined($Row->{"DATABASE ID"})) {
2206                                    #Getting the compound row
2207                                    my $Compartment = "c";
2208                                    if (defined($Row->{COMPARTMENT}->[0])) {
2209                                            $Compartment = $Row->{COMPARTMENT}->[0];
2210                                    }
2211                                    my $RxnRow = $rxnclasstable->get_row_by_key($Row->{"DATABASE ID"}->[0],"REACTION",1);
2212                                    my $Max = 0;
2213                                    my $Min = 0;
2214                                    my $Class = "Unknown";
2215                                    if (defined($DeadReactions) && defined($DeadReactions->{$Row->{"DATABASE ID"}->[0]})) {
2216                                            $Class = "Dead";
2217                                    } elsif (defined($Row->{"Min FLUX"}) && defined($Row->{"Max FLUX"})) {
2218                                            $Max = $Row->{"Max FLUX"}->[0];
2219                                            $Min = $Row->{"Min FLUX"}->[0];
2220                                            if ($Row->{"Min FLUX"}->[0] > 0.00000001) {
2221                                                    $Class = "Positive";
2222                                            } elsif ($Row->{"Max FLUX"}->[0] < -0.00000001) {
2223                                                    $Class = "Negative";
2224                                            } elsif ($Row->{"Min FLUX"}->[0] < -0.00000001) {
2225                                                    if ($Row->{"Max FLUX"}->[0] > 0.00000001) {
2226                                                            $Class = "Variable";
2227                                                    } else {
2228                                                            $Max = 0;
2229                                                            $Class = "Negative variable";
2230                                                    }
2231                                            } elsif ($Row->{"Max FLUX"}->[0] > 0.00000001) {
2232                                                    $Min = 0;
2233                                                    $Class = "Positive variable";
2234                                            } else {
2235                                                    $Min = 0;
2236                                                    $Max = 0;
2237                                                    $Class = "Blocked";
2238                                            }
2239                                    }
2240                                    my $index = 0;
2241                                    if (defined($RxnRow->{MEDIA})) {
2242                                            for (my $i=0; $i < @{$RxnRow->{MEDIA}};$i++) {
2243                                                    $index++;
2244                                                    if ($RxnRow->{MEDIA}->[$i] eq $Media) {
2245                                                            $index = $i;
2246                                                            last;
2247                                                    }
2248                                            }
2249                                    }
2250                                    $RxnRow->{MIN}->[$index] = $Min;
2251                                    $RxnRow->{MAX}->[$index] = $Max;
2252                                    $RxnRow->{CLASS}->[$index] = $Class;
2253                                    $RxnRow->{MEDIA}->[$index] = $Media;
2254                            }
2255                    }
2256                    if (!defined($SaveChanges) || $SaveChanges == 1) {
2257                            $rxnclasstable->save();
2258                    }
2259            }
2260            return ($rxnclasstable,$cpdclasstable);
2261    }
2262    
2263    =head3 RunAllStudiesWithDataFast
2264    Definition:
2265            (integer::false positives,integer::false negatives,integer::correct negatives,integer::correct positives,string::error vector,string heading vector) = FIGMODELmodel->RunAllStudiesWithDataFast(string::experiment,0/1::print result);
2266    Description:
2267            Simulates every experimental condition currently available for the model.
2268    =cut
2269    
2270    sub RunAllStudiesWithDataFast {
2271            my ($self,$Experiment,$PrintResults) = @_;
2272    
2273            #Printing lp and key file for model
2274            if (!-e $self->directory()."FBA-".$self->id().$self->selected_version().".lp") {
2275                    $self->PrintModelLPFile();
2276            }
2277            my $UniqueFilename = $self->figmodel()->filename();
2278    
2279            #Determing the simulations that need to be run
2280            my $ExperimentalDataTable = $self->figmodel()->GetExperimentalDataTable($self->genome(),$Experiment);
2281            #Creating the table of jobs to submit
2282            my $JobArray = $self->GetSimulationJobTable($ExperimentalDataTable,$Experiment,$UniqueFilename);
2283            #Printing the job file
2284            if (!-d $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/") {
2285                    system("mkdir ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/");
2286            }
2287            $JobArray->save();
2288    
2289            #Running simulations
2290            system($self->config("mfalite executable")->[0]." ".$self->config("Reaction database directory")->[0]."masterfiles/MediaTable.txt ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Jobfile.txt ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt");
2291            print "simulation complete for ".$self->id().$self->selected_version()."\n";
2292            #Parsing the results
2293            my $Results = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt",";","\\|",0,undef);
2294            if (!defined($Results)) {
2295                    $self->figmodel()->error_message("FIGMODELmodel:RunAllStudiesWithDataFast:Could not find simulation results: ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt");
2296                    return undef;
2297            }
2298            my ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector,$SimulationResults) = $self->EvaluateSimulationResults($Results,$ExperimentalDataTable);
2299            #Printing results to file
2300            $self->figmodel()->database()->save_table($SimulationResults,undef,undef,undef,"False negatives\tFalse positives\tCorrect negatives\tCorrect positives\n".$FalseNegatives."\t".$FalsePostives."\t".$CorrectNegatives."\t".$CorrectPositives."\n");
2301            $self->figmodel()->clearing_output($UniqueFilename);
2302    
2303            return ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector);
2304    }
2305    
2306    =head3 GetSimulationJobTable
2307    Definition:
2308            my $JobTable = $model->GetSimulationJobTable($Experiment,$PrintResults,$Version);
2309    Description:
2310    =cut
2311    
2312    sub GetSimulationJobTable {
2313            my ($self,$SimulationTable,$Experiment,$Folder) = @_;
2314    
2315            #Determing the simulations that need to be run
2316            if (!defined($SimulationTable)) {
2317                    $SimulationTable = $self->figmodel()->GetExperimentalDataTable($self->genome(),$Experiment);
2318                    if (!defined($SimulationTable)) {
2319                            return undef;
2320                    }
2321            }
2322    
2323            #Creating the job table
2324            my $JobTable = $self->figmodel()->CreateJobTable($Folder);
2325            for (my $i=0; $i < $SimulationTable->size(); $i++) {
2326                    if ($SimulationTable->get_row($i)->{"Heading"}->[0] =~ m/Gene\sKO/) {
2327                            my $Row = $JobTable->get_row_by_key("Gene KO","LABEL",1);
2328                            $JobTable->add_data($Row,"MEDIA",$SimulationTable->get_row($i)->{"Media"}->[0],1);
2329                    } elsif ($SimulationTable->get_row($i)->{"Heading"}->[0] =~ m/Media\sgrowth/) {
2330                            my $Row = $JobTable->get_row_by_key("Growth phenotype","LABEL",1);
2331                            $JobTable->add_data($Row,"MEDIA",$SimulationTable->get_row($i)->{"Media"}->[0],1);
2332                    } elsif ($SimulationTable->get_row($i)->{"Heading"}->[0] =~ m/Interval\sKO/) {
2333                            my $Row = $JobTable->get_row_by_key($SimulationTable->get_row($i)->{"Heading"}->[0],"LABEL",1);
2334                            $JobTable->add_data($Row,"MEDIA",$SimulationTable->get_row($i)->{"Media"}->[0],1);
2335                            $JobTable->add_data($Row,"GENE KO",$SimulationTable->get_row($i)->{"Experiment type"}->[0],1);
2336                    }
2337            }
2338    
2339            #Filling in model specific elements of the job table
2340            for (my $i=0; $i < $JobTable->size(); $i++) {
2341                    if ($JobTable->get_row($i)->{"LABEL"}->[0] =~ m/Gene\sKO/) {
2342                            $JobTable->get_row($i)->{"RUNTYPE"}->[0] = "SINGLEKO";
2343                            $JobTable->get_row($i)->{"SAVE NONESSENTIALS"}->[0] = 1;
2344                    } else {
2345                            $JobTable->get_row($i)->{"RUNTYPE"}->[0] = "GROWTH";
2346                            $JobTable->get_row($i)->{"SAVE NONESSENTIALS"}->[0] = 0;
2347                    }
2348                    $JobTable->get_row($i)->{"LP FILE"}->[0] = $self->directory()."FBA-".$self->id().$self->selected_version();
2349                    $JobTable->get_row($i)->{"MODEL"}->[0] = $self->directory().$self->id().$self->selected_version().".txt";
2350                    $JobTable->get_row($i)->{"SAVE FLUXES"}->[0] = 0;
2351            }
2352    
2353            return $JobTable;
2354    }
2355    
2356    =head3 EvaluateSimulationResults
2357    Definition:
2358            (integer::false positives,integer::false negatives,integer::correct negatives,integer::correct positives,string::error vector,string heading vector,FIGMODELtable::simulation results) = FIGMODELmodel->EvaluateSimulationResults(FIGMODELtable::raw simulation results,FIGMODELtable::experimental data);
2359    Description:
2360            Compares simulation results with experimental data to produce a table indicating where predictions are incorrect.
2361    =cut
2362    
2363    sub EvaluateSimulationResults {
2364            my ($self,$Results,$ExperimentalDataTable) = @_;
2365    
2366            #Comparing experimental results with simulation results
2367            my $SimulationResults = FIGMODELTable->new(["Run result","Experiment type","Media","Experiment ID","Reactions knocked out"],$self->directory()."SimulationOutput".$self->id().$self->selected_version().".txt",["Experiment ID","Media"],"\t",",",undef);
2368            my $FalsePostives = 0;
2369            my $FalseNegatives = 0;
2370            my $CorrectNegatives = 0;
2371            my $CorrectPositives = 0;
2372            my @Errorvector;
2373            my @HeadingVector;
2374            my $ReactionKOWithGeneHash;
2375            for (my $i=0; $i < $Results->size(); $i++) {
2376                    if ($Results->get_row($i)->{"LABEL"}->[0] eq "Gene KO") {
2377                            if (defined($Results->get_row($i)->{"REACTION KO WITH GENES"})) {
2378                                    for (my $j=0; $j < @{$Results->get_row($i)->{"REACTION KO WITH GENES"}}; $j++) {
2379                                            my @Temp = split(/:/,$Results->get_row($i)->{"REACTION KO WITH GENES"}->[$j]);
2380                                            if (defined($Temp[1]) && length($Temp[1]) > 0) {
2381                                                    $ReactionKOWithGeneHash->{$Temp[0]} = $Temp[1];
2382                                            }
2383                                    }
2384                            }
2385                            if ($Results->get_row($i)->{"OBJECTIVE"}->[0] == 0) {
2386                                    for (my $j=0; $j < @{$Results->get_row($i)->{"NONESSENTIALGENES"}}; $j++) {
2387                                            my $Row = $ExperimentalDataTable->get_row_by_key("Gene KO:".$Results->get_row($i)->{"MEDIA"}->[0].":".$Results->get_row($i)->{"NONESSENTIALGENES"}->[$j],"Heading");
2388                                            if (defined($Row)) {
2389                                                    my $KOReactions = "none";
2390                                                    if (defined($ReactionKOWithGeneHash->{$Results->get_row($i)->{"NONESSENTIALGENES"}->[$j]})) {
2391                                                            $KOReactions = $ReactionKOWithGeneHash->{$Results->get_row($i)->{"NONESSENTIALGENES"}->[$j]};
2392                                                    }
2393                                                    push(@HeadingVector,$Row->{"Heading"}->[0].":".$KOReactions);
2394                                                    my $Status = "Unknown";
2395                                                    if ($Row->{"Growth"}->[0] > 0) {
2396                                                            $Status = "False negative";
2397                                                            $FalseNegatives++;
2398                                                            push(@Errorvector,3);
2399                                                    } else {
2400                                                            $Status = "False positive";
2401                                                            $FalsePostives++;
2402                                                            push(@Errorvector,2);
2403                                                    }
2404                                                    $SimulationResults->add_row({"Run result" => [$Status],"Experiment type" => ["Gene KO"],"Media" => [$Row->{"Media"}->[0]],"Experiment ID" => [$Row->{"Experiment ID"}->[0]],"Reactions knocked out" => [$KOReactions]});
2405                                            }
2406                                    }
2407                            } else {
2408                                    for (my $j=0; $j < @{$Results->get_row($i)->{"ESSENTIALGENES"}}; $j++) {
2409                                            #print $j."\t".$Results->get_row($i)->{"ESSENTIALGENES"}->[$j]."\n";
2410                                            my $Row = $ExperimentalDataTable->get_row_by_key("Gene KO:".$Results->get_row($i)->{"MEDIA"}->[0].":".$Results->get_row($i)->{"ESSENTIALGENES"}->[$j],"Heading");
2411                                            if (defined($Row)) {
2412                                                    my $KOReactions = "none";
2413                                                    if (defined($ReactionKOWithGeneHash->{$Results->get_row($i)->{"ESSENTIALGENES"}->[$j]})) {
2414                                                            $KOReactions = $ReactionKOWithGeneHash->{$Results->get_row($i)->{"ESSENTIALGENES"}->[$j]};
2415                                                    }
2416                                                    push(@HeadingVector,$Row->{"Heading"}->[0].":".$KOReactions);
2417                                                    my $Status = "Unknown";
2418                                                    if ($Row->{"Growth"}->[0] > 0) {
2419                                                            $Status = "False negative";
2420                                                            $FalseNegatives++;
2421                                                            push(@Errorvector,3);
2422                                                    } else {
2423                                                            $Status = "Correct negative";
2424                                                            $CorrectNegatives++;
2425                                                            push(@Errorvector,1);
2426                                                    }
2427                                                    $SimulationResults->add_row({"Run result" => [$Status],"Experiment type" => ["Gene KO"],"Media" => [$Row->{"Media"}->[0]],"Experiment ID" => [$Row->{"Experiment ID"}->[0]],"Reactions knocked out" => [$KOReactions]});
2428                                            }
2429                                    }
2430                                    for (my $j=0; $j < @{$Results->get_row($i)->{"NONESSENTIALGENES"}}; $j++) {
2431                                            my $Row = $ExperimentalDataTable->get_row_by_key("Gene KO:".$Results->get_row($i)->{"MEDIA"}->[0].":".$Results->get_row($i)->{"NONESSENTIALGENES"}->[$j],"Heading");
2432                                            if (defined($Row)) {
2433                                                    my $KOReactions = "none";
2434                                                    if (defined($ReactionKOWithGeneHash->{$Results->get_row($i)->{"NONESSENTIALGENES"}->[$j]})) {
2435                                                            $KOReactions = $ReactionKOWithGeneHash->{$Results->get_row($i)->{"NONESSENTIALGENES"}->[$j]};
2436                                                    }
2437                                                    push(@HeadingVector,$Row->{"Heading"}->[0].":".$KOReactions);
2438                                                    my $Status = "Unknown";
2439                                                    if ($Row->{"Growth"}->[0] > 0) {
2440                                                            $Status = "Correct positive";
2441                                                            $CorrectPositives++;
2442                                                            push(@Errorvector,0);
2443                                                    } else {
2444                                                            $Status = "False positive";
2445                                                            $FalsePostives++;
2446                                                            push(@Errorvector,2);
2447                                                    }
2448                                                    $SimulationResults->add_row({"Run result" => [$Status],"Experiment type" => ["Gene KO"],"Media" => [$Row->{"Media"}->[0]],"Experiment ID" => [$Row->{"Experiment ID"}->[0]],"Reactions knocked out" => [$KOReactions]});
2449                                            }
2450                                    }
2451                            }
2452                    } elsif ($Results->get_row($i)->{"LABEL"}->[0] eq "Growth phenotype") {
2453                            my $Row = $ExperimentalDataTable->get_row_by_key("Media growth:".$Results->get_row($i)->{"MEDIA"}->[0].":".$Results->get_row($i)->{"MEDIA"}->[0],"Heading");
2454                            if (defined($Row)) {
2455                                    push(@HeadingVector,$Row->{"Heading"}->[0].":none");
2456                                    my $Status = "Unknown";
2457                                    if ($Row->{"Growth"}->[0] > 0) {
2458                                            if ($Results->get_row($i)->{"OBJECTIVE"}->[0] > 0) {
2459                                                    $Status = "Correct positive";
2460                                                    $CorrectPositives++;
2461                                                    push(@Errorvector,0);
2462                                            } else {
2463                                                    $Status = "False negative";
2464                                                    $FalseNegatives++;
2465                                                    push(@Errorvector,3);
2466                                            }
2467                                    } else {
2468                                            if ($Results->get_row($i)->{"OBJECTIVE"}->[0] > 0) {
2469                                                    $Status = "False positive";
2470                                                    $FalsePostives++;
2471                                                    push(@Errorvector,2);
2472                                            } else {
2473                                                    $Status = "Correct negative";
2474                                                    $CorrectNegatives++;
2475                                                    push(@Errorvector,1);
2476                                            }
2477                                    }
2478                                    $SimulationResults->add_row({"Run result" => [$Status],"Experiment type" => ["Media growth"],"Media" => [$Row->{"Media"}->[0]],"Experiment ID" => [$Row->{"Media"}->[0]],"Reactions knocked out" => ["none"]});
2479                            }
2480                    } elsif ($Results->get_row($i)->{"LABEL"}->[0] =~ m/Interval\sKO/ && defined($Results->get_row($i)->{"KOGENES"}->[0])) {
2481                            my $Row = $ExperimentalDataTable->get_row_by_key($Results->get_row($i)->{"LABEL"}->[0],"Heading");
2482                            if (defined($Row)) {
2483                                    my $Status = "Unknown";
2484                                    if ($Row->{"Growth"}->[0] > 0) {
2485                                            if ($Results->get_row($i)->{"OBJECTIVE"}->[0] > 0) {
2486                                                    $Status = "Correct positive";
2487                                                    $CorrectPositives++;
2488                                                    push(@Errorvector,0);
2489                                            } else {
2490                                                    $Status = "False negative";
2491                                                    $FalseNegatives++;
2492                                                    push(@Errorvector,3);
2493                                            }
2494                                    } else {
2495                                            if ($Results->get_row($i)->{"OBJECTIVE"}->[0] > 0) {
2496                                                    $Status = "False positive";
2497                                                    $FalsePostives++;
2498                                                    push(@Errorvector,2);
2499                                            } else {
2500                                                    $Status = "Correct negative";
2501                                                    $CorrectNegatives++;
2502                                                    push(@Errorvector,1);
2503                                            }
2504                                    }
2505                                    $SimulationResults->add_row({"Run result" => [$Status],"Experiment type" => ["Interval KO"],"Media" => [$Row->{"Media"}->[0]],"Experiment ID" => [$Row->{"Experiment ID"}->[0]],"Reactions knocked out" => ["none"]});
2506                            }
2507                    }
2508            }
2509    
2510            return ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,join(";",@Errorvector),join(";",@HeadingVector),$SimulationResults);
2511    }
2512    
2513    =head3 InspectSolution
2514    Definition:
2515            $model->InspectSolution(string::gene knocked out,string::media condition,[string]::list of reactions);
2516    Description:
2517    =cut
2518    
2519    sub InspectSolution {
2520            my ($self,$GeneKO,$Media,$ReactionList) = @_;
2521    
2522            #Getting a directory for the results
2523            my $UniqueFilename = $self->figmodel()->filename();
2524            system("mkdir ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/");
2525            my $TempVersion = "V".$UniqueFilename;
2526    
2527            #Setting gene ko to none if no genes are to be knocked out
2528            if ($GeneKO !~ m/^peg\./) {
2529                    $GeneKO = "none";
2530            }
2531    
2532            #Implementing the input solution in the test model
2533            my $ReactionArray;
2534            my $DirectionArray;
2535            my %SolutionHash;
2536            for (my $k=0; $k < @{$ReactionList}; $k++) {
2537                    if ($ReactionList->[$k] =~ m/(.+)(rxn\d\d\d\d\d)/) {
2538                            my $Reaction = $2;
2539                            my $Sign = $1;
2540                            if (defined($SolutionHash{$Reaction})) {
2541                                    $SolutionHash{$Reaction} = "<=>";
2542                            } elsif ($Sign eq "-") {
2543                                    $SolutionHash{$Reaction} = "<=";
2544                            } elsif ($Sign eq "+") {
2545                                    $SolutionHash{$Reaction} = "=>";
2546                            } else {
2547                                    $SolutionHash{$Reaction} = $Sign;
2548                            }
2549                    }
2550            }
2551            my @TempList = keys(%SolutionHash);
2552            for (my $k=0; $k < @TempList; $k++) {
2553                    push(@{$ReactionArray},$TempList[$k]);
2554                    push(@{$DirectionArray},$SolutionHash{$TempList[$k]});
2555            }
2556    
2557            print "Integrating solution!\n";
2558            $self->figmodel()->IntegrateGrowMatchSolution($self->id().$self->selected_version(),$self->directory().$self->id().$TempVersion.".txt",$ReactionArray,$DirectionArray,"SolutionInspection",1,1);
2559    
2560            #Printing lp and key file for model
2561            $self->PrintModelLPFile();
2562    
2563            #Running FBA on the test model
2564            my $JobTable = $self->figmodel()->CreateJobTable($UniqueFilename);
2565            $JobTable->add_row({"LABEL" => ["TEST"],"RUNTYPE" => ["GROWTH"],"LP FILE" => [$self->directory()."FBA-".$self->id().$TempVersion],"MODEL" => [$self->directory().$self->id().$TempVersion.".txt"],"MEDIA" => [$Media],"REACTION KO" => ["none|".join("|",@{$ReactionList})],"GENE KO" => [$GeneKO],"SAVE FLUXES" => [0],"SAVE NONESSENTIALS" => [0]});
2566            $JobTable->save();
2567    
2568            #Running simulations
2569            system($self->config("mfalite executable")->[0]." ".$self->config("Reaction database directory")->[0]."masterfiles/MediaTable.txt ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Jobfile.txt ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt");
2570    
2571            #Parsing the results
2572            my $Results = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt",";","\\|",0,undef);
2573            if (!defined($Results)) {
2574                    $self->figmodel()->error_message("FIGMODELmodel:InspectSolution:Could not load problem report ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt");
2575                    return undef;
2576            }
2577    
2578            #Making sure that the model grew with all reactions present
2579            my $Found = 0;
2580            for (my $i=0; $i < $Results->size(); $i++) {
2581                    if (defined($Results->get_row($i)->{"KOGENES"}->[0]) && defined($Results->get_row($i)->{"KOREACTIONS"}->[0]) && $Results->get_row($i)->{"KOREACTIONS"}->[0] eq "none" && $Results->get_row($i)->{"KOGENES"}->[0] eq $GeneKO && $Results->get_row($i)->{"OBJECTIVE"}->[0] > 0.00001) {
2582                            $Found = 1;
2583                    }
2584            }
2585            if ($Found == 0) {
2586                    print "Solution no longer valid\n";
2587                    return undef;
2588            }
2589    
2590            #Making sure all of the reactions added are still necessary
2591            my $FinalReactionList;
2592            for (my $k=0; $k < $Results->size(); $k++) {
2593                    if (defined($Results->get_row($k)->{"KOGENES"}->[0]) && $Results->get_row($k)->{"KOGENES"}->[0] eq $GeneKO) {
2594                            if (defined($Results->get_row($k)->{"KOREACTIONS"}->[0]) && $Results->get_row($k)->{"KOREACTIONS"}->[0] =~ m/rxn\d\d\d\d\d/ && $Results->get_row($k)->{"OBJECTIVE"}->[0] < 0.000001) {
2595                                    push(@{$FinalReactionList},$Results->get_row($k)->{"KOREACTIONS"}->[0]);
2596                            }
2597                    }
2598            }
2599    
2600            #Deleting extra files created
2601            unlink($self->directory()."FBA-".$self->id().$TempVersion.".lp");
2602            unlink($self->directory()."FBA-".$self->id().$TempVersion.".key");
2603            unlink($self->directory().$self->id().$TempVersion.".txt");
2604    
2605            #Deleting the test model and the MFA folder
2606            $self->figmodel()->clearing_output($UniqueFilename);
2607    
2608            return $FinalReactionList;
2609    }
2610    
2611    =head3 GapFillingAlgorithm
2612    
2613    Definition:
2614            FIGMODELmodel->GapFillingAlgorithm();
2615    
2616    Description:
2617            This is a wrapper for running the gap filling algorithm on any model in the database.
2618            The algorithm performs a gap filling for any false negative prediction of the avialable experimental data.
2619            This function is threaded to improve efficiency: one thread does nothing but using the MFAToolkit to fill gaps for every false negative prediction.
2620            The other thread reads in the gap filling solutions, builds a test model for each solution, and runs the test model against all available experimental data.
2621            This function prints two important output files in the Model directory:
2622            1.) GapFillingOutput.txt: this is a summary of the results of the gap filling analysis
2623            2.) GapFillingErrorMatrix.txt: this lists the correct and incorrect predictions for each gapfilling solution implemented in a test model.
2624    =cut
2625    
2626    sub GapFillingAlgorithm {
2627            my ($self) = @_;
2628    
2629            #First the input model version and model filename should be simulated and the false negatives identified
2630            my ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector) = $self->RunAllStudiesWithDataFast("All");
2631    
2632            #Getting the filename
2633            my $UniqueFilename = $self->figmodel()->filename();
2634    
2635            #Printing the original performance vector
2636            $self->figmodel()->database()->print_array_to_file($self->directory().$self->id().$self->selected_version()."-OPEM".".txt",[$HeadingVector,$Errorvector]);
2637    
2638            my $PreviousGapFilling;
2639            if (-e $self->directory().$self->id().$self->selected_version()."-GFS.txt") {
2640                    #Backing up the old solution file
2641                    system("cp ".$self->directory().$self->id().$self->selected_version()."-GFS.txt ".$self->directory().$self->id().$self->selected_version()."-OldGFS.txt");
2642                    unlink($self->directory().$self->id().$self->selected_version()."-GFS.txt");
2643            }
2644            if (-e $self->directory().$self->id().$self->selected_version()."-OldGFS.txt") {
2645                    #Reading in the solution file from the previous gap filling if it exists
2646                    $PreviousGapFilling = $self->figmodel()->database()->load_table($self->directory().$self->id().$self->selected_version()."-OldGFS.txt",";",",",0,["Experiment"]);
2647            }
2648    
2649            #Now we use the simulation output to make the gap filling run data
2650            my @Errors = split(/;/,$Errorvector);
2651            my @Headings = split(/;/,$HeadingVector);
2652            my $GapFillingRunSpecs = "";
2653            my $Count = 0;
2654            my $RescuedPreviousResults;
2655            my $RunCount = 0;
2656            my $SolutionExistedCount = 0;
2657            my $AcceptedSolutions = 0;
2658            my $RejectedSolutions = 0;
2659            my $NoExistingSolutions = 0;
2660            for (my $i=0; $i < @Errors; $i++) {
2661                    if ($Errors[$i] == 3) {
2662                            my @HeadingDataArray = split(/:/,$Headings[$i]);
2663                            if ($HeadingDataArray[2] !~ m/^peg\./ || $HeadingDataArray[3] ne "none") {
2664                                    my $SolutionFound = 0;
2665                                    if (defined($PreviousGapFilling) && defined($PreviousGapFilling->get_row_by_key($HeadingDataArray[2],"Experiment"))) {
2666                                            my @Rows = $PreviousGapFilling->get_rows_by_key($HeadingDataArray[2],"Experiment");
2667                                            for (my $j=0; $j < @Rows; $j++) {
2668                                                    if ($HeadingDataArray[2] =~ m/^peg\./) {
2669                                                            my $ReactionList = $self->InspectSolution($HeadingDataArray[2],$HeadingDataArray[1],$Rows[$j]->{"Solution reactions"});
2670                                                            if (defined($ReactionList)) {
2671                                                                    print join(",",@{$Rows[$j]->{"Solution reactions"}})."\t".join(",",@{$ReactionList})."\n";
2672                                                                    $SolutionFound++;
2673                                                                    push(@{$RescuedPreviousResults},$Rows[$j]->{"Experiment"}->[0].";".$Rows[$j]->{"Solution index"}->[0].";".$Rows[$j]->{"Solution cost"}->[0].";".join(",",@{$ReactionList}));
2674                                                                    $AcceptedSolutions++;
2675                                                            } else {
2676                                                                    $RejectedSolutions++;
2677                                                            }
2678                                                    } else {
2679                                                            my $ReactionList = $self->InspectSolution($HeadingDataArray[2],$HeadingDataArray[1],$Rows[$j]->{"Solution reactions"});
2680                                                            if (defined($ReactionList)) {
2681                                                                    print join(",",@{$Rows[$j]->{"Solution reactions"}})."\t".join(",",@{$ReactionList})."\n";
2682                                                                    $SolutionFound++;
2683                                                                    push(@{$RescuedPreviousResults},$Rows[$j]->{"Experiment"}->[0].";".$Rows[$j]->{"Solution index"}->[0].";".$Rows[$j]->{"Solution cost"}->[0].";".join(",",@{$ReactionList}));
2684                                                                    $AcceptedSolutions++;
2685                                                            } else {
2686                                                                    $RejectedSolutions++;
2687                                                            }
2688                                                    }
2689                                            }
2690                                    } else {
2691                                            $NoExistingSolutions++;
2692                                    }
2693                                    if ($SolutionFound == 0) {
2694                                            $RunCount++;
2695                                            if (length($GapFillingRunSpecs) > 0) {
2696                                                    $GapFillingRunSpecs .= ";";
2697                                            }
2698                                            $GapFillingRunSpecs .= $HeadingDataArray[2].":".$HeadingDataArray[1].".txt:".$HeadingDataArray[3];
2699                                    } else {
2700                                            $SolutionExistedCount++;
2701                                    }
2702                            }
2703                            $Count++;
2704                    }
2705            }
2706    
2707            #Updating the growmatch progress table
2708            my $Row = $self->figmodel()->database()->get_row_by_key("GROWMATCH TABLE",$self->genome(),"ORGANISM",1);
2709            $Row->{"INITIAL FP"}->[0] = $FalsePostives;
2710            $Row->{"INITIAL FN"}->[0] = $FalseNegatives;
2711            $Row->{"GF TIMING"}->[0] = time()."-";
2712            $Row->{"FN WITH SOL"}->[0] = $FalseNegatives-$NoExistingSolutions;
2713            $Row->{"FN WITH ACCEPTED SOL"}->[0] = $SolutionExistedCount;
2714            $Row->{"TOTAL ACCEPTED GF SOL"}->[0] = $AcceptedSolutions;
2715            $Row->{"TOTAL REJECTED GF SOL"}->[0] = $RejectedSolutions;
2716            $Row->{"FN WITH NO SOL"}->[0] = $NoExistingSolutions+$RejectedSolutions;
2717            $self->figmodel()->database()->update_row("GROWMATCH TABLE",$Row,"ORGANISM");
2718    
2719            #Running the gap filling once to correct all false negative errors
2720            my $SolutionsFound = 0;
2721            my $GapFillingArray;
2722            push(@{$GapFillingArray},split(/;/,$GapFillingRunSpecs));
2723            $self->datagapfill($GapFillingArray);
2724    
2725            if (defined($RescuedPreviousResults) && @{$RescuedPreviousResults} > 0) {
2726                    #Printing previous solutions to GFS file
2727                    $self->figmodel()->database()->print_array_to_file($self->directory().$self->id().$self->selected_version()."-GFS.txt",$RescuedPreviousResults,1);
2728                    $SolutionsFound = 1;
2729            }
2730    
2731            #Recording the finishing of the gapfilling
2732            $Row = $self->figmodel()->database()->get_row_by_key("GROWMATCH TABLE",$self->genome(),"ORGANISM",1);
2733            $Row->{"GF TIMING"}->[0] .= time();
2734            $self->figmodel()->database()->update_row("GROWMATCH TABLE",$Row,"ORGANISM");
2735    
2736            if ($SolutionsFound == 1) {
2737                    #Scheduling solution testing
2738                    $self->figmodel()->add_job_to_queue("testsolutions?".$self->id().$self->selected_version()."?-1?GF","QSUB","fast","master","BACK");
2739            } else {
2740                    $self->figmodel()->error_message("No false negative predictions found. Data gap filling not necessary!");
2741            }
2742    
2743            return $self->success();
2744    }
2745    
2746    =head3 BuildSpecificBiomassReaction
2747    Definition:
2748            FIGMODELmodel->BuildSpecificBiomassReaction();
2749    Description:
2750    =cut
2751    sub BuildSpecificBiomassReaction {
2752            my ($self) = @_;
2753    
2754            my $biomassrxn;
2755            my $OrganismID = $self->genome();
2756            #Checking for a biomass override
2757            if (defined($self->config("biomass reaction override")->{$OrganismID})) {
2758                    $biomassrxn = $self->config("biomass reaction override")->{$OrganismID};
2759                    print "Overriding biomass template and selecting ".$biomassrxn." for ".$OrganismID.".\n";
2760            } else {#Creating biomass reaction from the template
2761                    #Getting the genome stats
2762                    my $genomestats = $self->figmodel()->get_genome_stats($self->genome());
2763                    my $Class = $genomestats->{CLASS}->[0];
2764                    my $Name = $genomestats->{NAME}->[0];
2765    
2766                    #Checking for phoenix variants
2767                    my $PhoenixVariantTable = $self->figmodel()->database()->GetDBTable("Phoenix variants table");
2768                    my $Phoenix = 0;
2769                    my @Rows = $PhoenixVariantTable->get_rows_by_key($OrganismID,"GENOME");
2770                    my $VariantHash;
2771                    for (my $i=0; $i < @Rows; $i++) {
2772                            $Phoenix = 1;
2773                            if (defined($Rows[$i]->{"SUBSYSTEM"}) && defined($Rows[$i]->{"VARIANT"})) {
2774                                    $VariantHash->{$Rows[$i]->{"SUBSYSTEM"}->[0]} = $Rows[$i]->{"VARIANT"}->[0];
2775                            }
2776                    }
2777    
2778                    #Collecting genome data
2779                    my $RoleHash;
2780                    my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome());
2781                    for (my $i=0; $i < $FeatureTable->size(); $i++) {
2782                            for (my $j=0; $j < @{$FeatureTable->get_row($i)->{"ROLES"}}; $j++) {
2783                                    $RoleHash->{$FeatureTable->get_row($i)->{"ROLES"}->[$j]} = 1;
2784                                    my $Subsystems = $self->figmodel()->subsystems_of_role($FeatureTable->get_row($i)->{"ROLES"}->[$j]);
2785                                    if (defined($Subsystems)) {
2786                                            for (my $k=0; $k < @{$Subsystems}; $k++) {
2787                                                    if ($Phoenix == 0) {
2788                                                            $VariantHash->{$Subsystems->[$k]} = 1;
2789                                                    }
2790                                            }
2791                                    }
2792                            }
2793                    }
2794    
2795                    #Scanning through the template item by item and determinine which biomass components should be added
2796                    my $ComponentTypes;
2797                    my $EquationData;
2798                    my $BiomassReactionTemplateTable = $self->figmodel()->database()->GetDBTable("BIOMASS TEMPLATE");
2799                    for (my $i=0; $i < $BiomassReactionTemplateTable->size(); $i++) {
2800                            my $Row = $BiomassReactionTemplateTable->get_row($i);
2801                            if (defined($Row->{"INCLUSION CRITERIA"}->[0]) && $Row->{"INCLUSION CRITERIA"}->[0] eq "UNIVERSAL") {
2802                                    push(@{$EquationData},$Row);
2803                                    $ComponentTypes->{$Row->{"CLASS"}->[0]}->{$Row->{"ID"}->[0]} = 1;
2804                            } else {
2805                                    my $Criteria = $Row->{"INCLUSION CRITERIA"}->[0];
2806                                    my $End = 0;
2807                                    while ($End == 0) {
2808                                            if ($Criteria =~ m/^(.+)(AND)\{([^{^}]+)\}(.+)$/ || $Criteria =~ m/^(AND)\{([^{^}]+)\}$/ || $Criteria =~ m/^(.+)(OR)\{([^{^}]+)\}(.+)$/ || $Criteria =~ m/^(OR)\{([^{^}]+)\}$/) {
2809                                                    print $Criteria."\n";
2810                                                    my $Start = "";
2811                                                    my $End = "";
2812                                                    my $Condition = $1;
2813                                                    my $Data = $2;
2814                                                    if ($1 ne "AND" && $1 ne "OR") {
2815                                                            $Start = $1;
2816                                                            $End = $4;
2817                                                            $Condition = $2;
2818                                                            $Data = $3;
2819                                                    }
2820                                                    my $Result = "YES";
2821                                                    if ($Condition eq "OR") {
2822                                                            $Result = "NO";
2823                                                    }
2824                                                    my @Array = split(/\|/,$Data);
2825                                                    for (my $j=0; $j < @Array; $j++) {
2826                                                            if ($Array[$j] eq "YES" && $Condition eq "OR") {
2827                                                                    $Result = "YES";
2828                                                                    last;
2829                                                            } elsif ($Array[$j] eq "NO" && $Condition eq "AND") {
2830                                                                    $Result = "NO";
2831                                                                    last;
2832                                                            } elsif ($Array[$j] =~ m/^COMPOUND:(.+)/) {
2833                                                                    my $Match = 0;
2834                                                                    for (my $k=0; $k < @{$EquationData}; $k++) {
2835                                                                            if ($EquationData->[$k]->{"ID"}->[0] eq $1) {
2836                                                                                    $Match = 1;
2837                                                                                    last;
2838                                                                            }
2839                                                                    }
2840                                                                    if ($Match == 1 && $Condition eq "OR") {
2841                                                                            $Result = "YES";
2842                                                                            last;
2843                                                                    } elsif ($Match != 1 && $Condition eq "AND") {
2844                                                                            $Result = "NO";
2845                                                                            last;
2846                                                                    }
2847                                                            } elsif ($Array[$j] =~ m/^!COMPOUND:(.+)/) {
2848                                                                    my $Match = 0;
2849                                                                    for (my $k=0; $k < @{$EquationData}; $k++) {
2850                                                                            if ($EquationData->[$k]->{"ID"}->[0] eq $1) {
2851                                                                                    $Match = 1;
2852                                                                                    last;
2853                                                                            }
2854                                                                    }
2855                                                                    if ($Match != 1 && $Condition eq "OR") {
2856                                                                            $Result = "YES";
2857                                                                            last;
2858                                                                    } elsif ($Match == 1 && $Condition eq "AND") {
2859                                                                            $Result = "NO";
2860                                                                            last;
2861                                                                    }
2862                                                            } elsif ($Array[$j] =~ m/^NAME:(.+)/) {
2863                                                                    my $Comparison = $1;
2864                                                                    if ((!defined($Comparison) || !defined($Name) || $Name =~ m/$Comparison/) && $Condition eq "OR") {
2865                                                                            $Result = "YES";
2866                                                                            last;
2867                                                                    } elsif (defined($Comparison) && defined($Name) && $Name !~ m/$Comparison/ && $Condition eq "AND") {
2868                                                                            $Result = "NO";
2869                                                                            last;
2870                                                                    }
2871                                                            } elsif ($Array[$j] =~ m/^!NAME:(.+)/) {
2872                                                                    my $Comparison = $1;
2873                                                                    if ((!defined($Comparison) || !defined($Name) || $Name !~ m/$Comparison/) && $Condition eq "OR") {
2874                                                                            $Result = "YES";
2875                                                                            last;
2876                                                                    } elsif (defined($Comparison) && defined($Name) && $Name =~ m/$Comparison/ && $Condition eq "AND") {
2877                                                                            $Result = "NO";
2878                                                                            last;
2879                                                                    }
2880                                                            } elsif ($Array[$j] =~ m/^SUBSYSTEM:(.+)/) {
2881                                                                    my @SubsystemArray = split(/`/,$1);
2882                                                                    if (@SubsystemArray == 1) {
2883                                                                            if (defined($VariantHash->{$SubsystemArray[0]}) && $VariantHash->{$SubsystemArray[0]} ne -1 && $Condition eq "OR") {
2884                                                                                    $Result = "YES";
2885                                                                                    last;
2886                                                                            } elsif ((!defined($VariantHash->{$SubsystemArray[0]}) || $VariantHash->{$SubsystemArray[0]} eq -1) && $Condition eq "AND") {
2887                                                                                    $Result = "NO";
2888                                                                                    last;
2889                                                                            }
2890                                                                    } else {
2891                                                                            my $Match = 0;
2892                                                                            for (my $k=1; $k < @SubsystemArray; $k++) {
2893                                                                                    if (defined($VariantHash->{$SubsystemArray[0]}) && $VariantHash->{$SubsystemArray[0]} eq $SubsystemArray[$k]) {
2894                                                                                            $Match = 1;
2895                                                                                            last;
2896                                                                                    }
2897                                                                            }
2898                                                                            if ($Match == 1 && $Condition eq "OR") {
2899                                                                                    $Result = "YES";
2900                                                                                    last;
2901                                                                            } elsif ($Match != 1 && $Condition eq "AND") {
2902                                                                                    $Result = "NO";
2903                                                                                    last;
2904                                                                            }
2905                                                                    }
2906                                                            } elsif ($Array[$j] =~ m/^!SUBSYSTEM:(.+)/) {
2907                                                                    my @SubsystemArray = split(/`/,$1);
2908                                                                    if (@SubsystemArray == 1) {
2909                                                                            if ((!defined($VariantHash->{$SubsystemArray[0]}) || $VariantHash->{$SubsystemArray[0]} eq -1) && $Condition eq "OR") {
2910                                                                                    $Result = "YES";
2911                                                                                    last;
2912                                                                            } elsif (defined($VariantHash->{$SubsystemArray[0]}) && $VariantHash->{$SubsystemArray[0]} ne -1 && $Condition eq "AND") {
2913                                                                                    $Result = "NO";
2914                                                                                    last;
2915                                                                            }
2916                                                                    } else {
2917                                                                            my $Match = 0;
2918                                                                            for (my $k=1; $k < @SubsystemArray; $k++) {
2919                                                                                    if (defined($VariantHash->{$SubsystemArray[0]}) && $VariantHash->{$SubsystemArray[0]} eq $SubsystemArray[$k]) {
2920                                                                                            $Match = 1;
2921                                                                                            last;
2922                                                                                    }
2923                                                                            }
2924                                                                            if ($Match != 1 && $Condition eq "OR") {
2925                                                                                    $Result = "YES";
2926                                                                                    last;
2927                                                                            } elsif ($Match == 1 && $Condition eq "AND") {
2928                                                                                    $Result = "NO";
2929                                                                                    last;
2930                                                                            }
2931                                                                    }
2932                                                            } elsif ($Array[$j] =~ m/^ROLE:(.+)/) {
2933                                                                    if (defined($RoleHash->{$1}) && $Condition eq "OR") {
2934                                                                            $Result = "YES";
2935                                                                            last;
2936                                                                    } elsif (!defined($RoleHash->{$1}) && $Condition eq "AND") {
2937                                                                            $Result = "NO";
2938                                                                            last;
2939                                                                    }
2940                                                            } elsif ($Array[$j] =~ m/^!ROLE:(.+)/) {
2941                                                                    if (!defined($RoleHash->{$1}) && $Condition eq "OR") {
2942                                                                            $Result = "YES";
2943                                                                            last;
2944                                                                    } elsif (defined($RoleHash->{$1}) && $Condition eq "AND") {
2945                                                                            $Result = "NO";
2946                                                                            last;
2947                                                                    }
2948                                                            } elsif ($Array[$j] =~ m/^CLASS:(.+)/) {
2949                                                                    if ($Class eq $1 && $Condition eq "OR") {
2950                                                                            $Result = "YES";
2951                                                                            last;
2952                                                                    } elsif ($Class ne $1 && $Condition eq "AND") {
2953                                                                            $Result = "NO";
2954                                                                            last;
2955                                                                    }
2956                                                            } elsif ($Array[$j] =~ m/^!CLASS:(.+)/) {
2957                                                                    if ($Class ne $1 && $Condition eq "OR") {
2958                                                                            $Result = "YES";
2959                                                                            last;
2960                                                                    } elsif ($Class eq $1 && $Condition eq "AND") {
2961                                                                            $Result = "NO";
2962                                                                            last;
2963                                                                    }
2964                                                            }
2965                                                    }
2966                                                    $Criteria = $Start.$Result.$End;
2967                                                    print $Criteria."\n";
2968                                            } else {
2969                                                    $End = 1;
2970                                                    last;
2971