[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.2, Fri Dec 11 21:37:02 2009 UTC revision 1.14, Thu May 13 05:04:26 2010 UTC
# Line 1  Line 1 
 package FIGMODELmodel;  
1  use strict;  use strict;
2    package FIGMODELmodel;
3    
4  =head1 FIGMODELmodel object  =head1 FIGMODELmodel object
5  =head2 Introduction  =head2 Introduction
# Line 14  Line 14 
14  =cut  =cut
15  sub new {  sub new {
16          my ($class,$figmodel,$id) = @_;          my ($class,$figmodel,$id) = @_;
17    
18          #Error checking first          #Error checking first
19          if (!defined($figmodel)) {          if (!defined($figmodel)) {
20                  print STDERR "FIGMODELmodel->new(undef,".$id."):figmodel must be defined to create a model object!\n";                  print STDERR "FIGMODELmodel->new(undef,".$id."):figmodel must be defined to create a model object!\n";
# Line 27  Line 28 
28          }          }
29    
30          #Checking that the id exists          #Checking that the id exists
31          my $tbl = $self->figmodel()->database()->GetDBTable("MODELS");          my $modelHandle = $self->figmodel()->database()->get_object_manager("model");
32          if (!defined($tbl)) {          if (!defined($modelHandle)) {
33                  $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,".$id."):could not load MODELS table. Check database!");                  $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,".$id."):could not load MODELS table. Check database!");
34                  return undef;                  return undef;
35          }          }
36    
37          #If the id is a number, we get the model row by index          #If the id is a number, we get the model row by index
         my $index = $id;  
38          if ($id =~ m/^\d+$/) {          if ($id =~ m/^\d+$/) {
39                  $self->{_data} = $tbl->get_row($id);                  my $objects = $modelHandle->get_objects();
40          } else {                  $self->{_data} = $objects->[$id];
41                  $self->{_data} = $tbl->get_row_by_key($id,"id");                  $self->figmodel()->{_models}->{$id} = $self;
42                  if (defined($self->{_data})) {          } else {
43                          $index = $tbl->row_index($self->{_data});                  my $objects = $modelHandle->get_objects({id => $id});
44                    if (!defined($objects->[0]) && $id =~ m/(.+)(V[^V]+)/) {
45                            $objects = $modelHandle->get_objects({id => $1});
46                            if (!defined($objects->[0])) {
47                                    $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,".$id."):could not find model ".$id." in database!");
48                                    return undef;
49                  }                  }
50                            $self->{_selectedversion} = $2;
51                    } elsif (!defined($objects->[0])) {
52                            $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,".$id."):could not find model ".$id." in database!");
53                            return undef;
54                    }
55                    $self->{_data} = $objects->[0];
56          }          }
57          if (!defined($self->{_data})) {          if (!defined($self->{_data})) {
58                  $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,".$id."):could not find specified id in database!");                  $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,".$id."):could not find specified id in database!");
59                  return undef;                  return undef;
60          }          }
         $self->{_index} = $index;  
61          $self->figmodel()->{_models}->{$self->id()} = $self;          $self->figmodel()->{_models}->{$self->id()} = $self;
         $self->figmodel()->{_models}->{$self->index()} = $self;  
62    
63          return $self;          return $self;
64  }  }
# Line 107  Line 116 
116  =cut  =cut
117  sub fig {  sub fig {
118          my ($self) = @_;          my ($self) = @_;
   
119          if (!defined($self->{_fig}) && $self->source() !~ /^MGRAST/) {          if (!defined($self->{_fig}) && $self->source() !~ /^MGRAST/) {
120                  if ($self->source() =~ /^RAST/) {                  if ($self->source() =~ /^RAST/) {
121                          $self->{"_fig"} = $self->figmodel()->fig();#$self->genome());                          $self->{"_fig"} = $self->figmodel()->fig();#$self->genome());
# Line 115  Line 123 
123                          $self->{"_fig"} = $self->figmodel()->fig();                          $self->{"_fig"} = $self->figmodel()->fig();
124                  }                  }
125          }          }
   
126          return $self->{"_fig"};          return $self->{"_fig"};
127  }  }
128    
129    =head3 aquireModelLock
130    
131    Definition:
132    
133            FIGMODELmodel->aquireModelLock();
134    
135    Description:
136    
137            Locks the database for alterations relating to the current model object
138    
139    =cut
140    sub aquireModelLock {
141            my ($self) = @_;
142            $self->figmodel()->database()->genericLock($self->id());
143    }
144    
145    =head3 releaseModelLock
146    
147    Definition:
148    
149            FIGMODELmodel->releaseModelLock();
150    
151    Description:
152    
153            Unlocks the database for alterations relating to the current model object
154    
155    =cut
156    sub releaseModelLock {
157            my ($self) = @_;
158            $self->figmodel()->database()->genericUnlock($self->id());
159    }
160    
161  =head3 mgdata  =head3 mgdata
162  Definition:  Definition:
163          FIGMODEL = FIGMODELmodel->mgdata();          FIGMODEL = FIGMODELmodel->mgdata();
# Line 127  Line 166 
166  =cut  =cut
167  sub mgdata {  sub mgdata {
168          my ($self) = @_;          my ($self) = @_;
   
169          if (!defined($self->{_mgdata}) && $self->source() =~ /^MGRAST/) {          if (!defined($self->{_mgdata}) && $self->source() =~ /^MGRAST/) {
170                  require MGRAST;                  require MGRAST;
171                  $self->{_mgdata} = $self->figmodel()->mgrast()->Job->get_objects( { 'genome_id' => $self->genome() } )                  $self->{_mgdata} = $self->figmodel()->mgrast()->Job->get_objects( { 'genome_id' => $self->genome() } )
172          }          }
   
173          return $self->{_mgdata};          return $self->{_mgdata};
174  }  }
175    
# Line 144  Line 181 
181  =cut  =cut
182  sub id {  sub id {
183          my ($self) = @_;          my ($self) = @_;
184          return $self->{_data}->{id}->[0];          return $self->{_data}->id();
185  }  }
186    
187  =head3 owner  =head3 owner
# Line 155  Line 192 
192  =cut  =cut
193  sub owner {  sub owner {
194          my ($self) = @_;          my ($self) = @_;
195          return $self->{_data}->{owner}->[0];          return $self->{_data}->owner();
 }  
   
 =head3 index  
 Definition:  
         string = FIGMODELmodel->index();  
 Description:  
         Returns model index  
 =cut  
 sub index {  
         my ($self) = @_;  
         return $self->{_index};  
196  }  }
197    
198  =head3 name  =head3 name
# Line 185  Line 211 
211                          if (defined($self->mgdata())) {                          if (defined($self->mgdata())) {
212                                  $self->{_name} = $self->mgdata()->genome_name;                                  $self->{_name} = $self->mgdata()->genome_name;
213                          }                          }
                 } elsif (defined($self->stats())) {  
                         $self->{_name} = $self->stats()->{'Organism name'}->[0];  
214                  } elsif ($source !~ /^RAST/) {                  } elsif ($source !~ /^RAST/) {
215                          $self->{_name} = $self->fig()->orgname_of_orgid($self->genome());                          $self->{_name} = $self->fig()->orgname_of_orgid($self->genome());
216                    } else {
217                            $self->{_name} = $self->figmodel()->get_genome_stats($self->genome())->{NAME}->[0];
218                  }                  }
219          }          }
220    
221          return $self->{_name};          return $self->{_name};
222  }  }
223    
 =head3 stats  
 Definition:  
         string = FIGMODELmodel->stats();  
 Description:  
         Returns model stats  
 =cut  
 sub stats {  
         my ($self) = @_;  
   
         if (!defined($self->{_stats})) {  
                 $self->{_stats} = $self->figmodel()->database()->GetDBTable("MODEL STATS")->get_row_by_key($self->id(),"Model ID");  
         }  
         return $self->{_stats};  
 }  
   
224  =head3 get_reaction_class  =head3 get_reaction_class
225  Definition:  Definition:
226          string = FIGMODELmodel->get_reaction_class(string::reaction ID);          string = FIGMODELmodel->get_reaction_class(string::reaction ID);
# Line 230  Line 241 
241                  my $ClassRow = $self->{_reaction_classes}->get_row_by_key($reaction,"REACTION");                  my $ClassRow = $self->{_reaction_classes}->get_row_by_key($reaction,"REACTION");
242                  if (defined($ClassRow) && defined($ClassRow->{CLASS})) {                  if (defined($ClassRow) && defined($ClassRow->{CLASS})) {
243                          my $class;                          my $class;
244                            my $min = $ClassRow->{MIN}->[0];
245                            my $max = $ClassRow->{MAX}->[0];
246                          if ($ClassRow->{CLASS}->[0] eq "Positive") {                          if ($ClassRow->{CLASS}->[0] eq "Positive") {
247                                  $class = "Essential =>";                                  $class = "Essential =>";
248                                    $class.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>";
249                          } elsif ($ClassRow->{CLASS}->[0] eq "Negative") {                          } elsif ($ClassRow->{CLASS}->[0] eq "Negative") {
250                                  $class = "Essential <=";                                  $class = "Essential <=";
251                                    $class.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>";
252                          } elsif ($ClassRow->{CLASS}->[0] eq "Positive variable") {                          } elsif ($ClassRow->{CLASS}->[0] eq "Positive variable") {
253                                  $class = "Active =>";                                  $class = "Active =>";
254                                    $class.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>";
255                          } elsif ($ClassRow->{CLASS}->[0] eq "Negative variable") {                          } elsif ($ClassRow->{CLASS}->[0] eq "Negative variable") {
256                                  $class = "Active <=";                                  $class = "Active <=";
257                                    $class.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>";
258                          } elsif ($ClassRow->{CLASS}->[0] eq "Variable") {                          } elsif ($ClassRow->{CLASS}->[0] eq "Variable") {
259                                  $class = "Active <=>";                                  $class = "Active <=>";
260                                    $class.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>";
261                          } elsif ($ClassRow->{CLASS}->[0] eq "Blocked") {                          } elsif ($ClassRow->{CLASS}->[0] eq "Blocked") {
262                                  $class = "Inactive";                                  $class = "Inactive";
263                          } elsif ($ClassRow->{CLASS}->[0] eq "Dead") {                          } elsif ($ClassRow->{CLASS}->[0] eq "Dead") {
264                                  $class = "Disconnected";                                  $class = "Disconnected";
265                          }                          }
266    
267                          if (!defined($nohtml) || $nohtml ne "1") {                          if (!defined($nohtml) || $nohtml ne "1") {
268                                  $class = "<span title=\"Flux:".$ClassRow->{MIN}->[0]." to ".$ClassRow->{MAX}->[0]."\">".$class."</span>";                                  $class = "<span title=\"Flux:".$min." to ".$max."\">".$class."</span>";
269                          }                          }
270    
271                          return $class;                          return $class;
272                  }                  }
273                  return undef;                  return undef;
# Line 268  Line 288 
288                                  $classstring .= "<br>";                                  $classstring .= "<br>";
289                          }                          }
290                          my $NewClass;                          my $NewClass;
291                            my $min = $ClassRow->{MIN}->[$i];
292                            my $max = $ClassRow->{MAX}->[$i];
293                          if ($ClassRow->{CLASS}->[$i] eq "Positive") {                          if ($ClassRow->{CLASS}->[$i] eq "Positive") {
294                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Essential =>";                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Essential =>";
295                                  if (!defined($nohtml) || $nohtml ne "1") {                                  $NewClass.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>";
                                         $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";  
                                 }  
296                          } elsif ($ClassRow->{CLASS}->[$i] eq "Negative") {                          } elsif ($ClassRow->{CLASS}->[$i] eq "Negative") {
297                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Essential <=";                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Essential <=";
298                                  if (!defined($nohtml) || $nohtml ne "1") {                                  $NewClass.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>";
                                         $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";  
                                 }  
299                          } elsif ($ClassRow->{CLASS}->[$i] eq "Positive variable") {                          } elsif ($ClassRow->{CLASS}->[$i] eq "Positive variable") {
300                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Active =>";                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Active =>";
301                                  if (!defined($nohtml) || $nohtml ne "1") {                                  $NewClass.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>";
                                         $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";  
                                 }  
302                          } elsif ($ClassRow->{CLASS}->[$i] eq "Negative variable") {                          } elsif ($ClassRow->{CLASS}->[$i] eq "Negative variable") {
303                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Active <=";                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Active <=";
304                                  if (!defined($nohtml) || $nohtml ne "1") {                                  $NewClass.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>";
                                         $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";  
                                 }  
305                          } elsif ($ClassRow->{CLASS}->[$i] eq "Variable") {                          } elsif ($ClassRow->{CLASS}->[$i] eq "Variable") {
306                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Active <=>";                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Active <=>";
307                                  if (!defined($nohtml) || $nohtml ne "1") {                                  $NewClass.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>";
                                         $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";  
                                 }  
308                          } elsif ($ClassRow->{CLASS}->[$i] eq "Blocked") {                          } elsif ($ClassRow->{CLASS}->[$i] eq "Blocked") {
309                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Inactive";                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Inactive";
                                 if (!defined($nohtml) || $nohtml ne "1") {  
                                         $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";  
                                 }  
310                          } elsif ($ClassRow->{CLASS}->[$i] eq "Dead") {                          } elsif ($ClassRow->{CLASS}->[$i] eq "Dead") {
311                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Disconnected";                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Disconnected";
                                 if (!defined($nohtml) || $nohtml ne "1") {  
                                         $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";  
312                                  }                                  }
313    
314                            if (!defined($nohtml) || $nohtml ne "1") {
315                                    $NewClass = "<span title=\"Flux:".$min." to ".$max."\">".$NewClass."</span>";
316                          }                          }
317                          $classstring .= $NewClass;                          $classstring .= $NewClass;
318                  }                  }
# Line 366  Line 376 
376          return $Data->{LOAD}->[0];          return $Data->{LOAD}->[0];
377  }  }
378    
379    =head3 load_model_table
380    
381    Definition: FIGMODELTable = FIGMODELmodel->load_model_table(string:table name,0/1:refresh the table));
382    
383    Description: Returns the table specified by the input filename. Table will be stored in a file in the model directory.
384    
385    =cut
386    sub load_model_table {
387            my ($self,$name,$refresh) = @_;
388            if (!defined($refresh)) {
389                    $refresh = 1;
390            }
391            if ($refresh == 1) {
392                    delete $self->{"_".$name};
393            }
394            if (!defined($self->{"_".$name})) {
395                    my $tbldef = $self->figmodel()->config("$name");
396                    if (!defined($tbldef)) {
397                            return undef;
398                    }
399                    my $filename = $self->directory().$name."-".$self->id().$self->selected_version().".tbl";
400                    $self->{"_".$name} = $self->figmodel()->database()->load_table($filename,"\t","|",$tbldef->{headingline}->[0],$tbldef->{hashcolumns});
401                    if (!defined($self->{"_".$name})) {
402                            if (defined($tbldef->{prefix})) {
403                                    $self->{"_".$name} = FIGMODELTable->new($tbldef->{columns},$filename,$tbldef->{hashcolumns},"\t","|",join(@{$tbldef->{prefix}},"\n"));
404                            } else {
405                                    $self->{"_".$name} = FIGMODELTable->new($tbldef->{columns},$filename,$tbldef->{hashcolumns},"\t","|");
406                            }
407                    }
408            }
409            return $self->{"_".$name};
410    }
411    
412  =head3 get_reaction_number  =head3 get_reaction_number
413  Definition:  Definition:
414          int = FIGMODELmodel->get_reaction_number();          int = FIGMODELmodel->get_reaction_number();
# Line 374  Line 417 
417  =cut  =cut
418  sub get_reaction_number {  sub get_reaction_number {
419          my ($self) = @_;          my ($self) = @_;
   
420          if (!defined($self->reaction_table())) {          if (!defined($self->reaction_table())) {
421                  return 0;                  return 0;
422          }          }
   
423          return $self->reaction_table()->size();          return $self->reaction_table()->size();
424  }  }
425    
# Line 393  Line 434 
434    
435          if (!defined($self->{_reaction_data})) {          if (!defined($self->{_reaction_data})) {
436                  $self->{_reaction_data} = $self->figmodel()->database()->GetDBModel($self->id());                  $self->{_reaction_data} = $self->figmodel()->database()->GetDBModel($self->id());
437                    my $classTbl = $self->reaction_class_table();
438                    for (my $i=0; $i < $classTbl->size(); $i++) {
439                            my $row = $classTbl->get_row($i);
440                            my $rxnRow = $self->{_reaction_data}->get_row_by_key($row->{"REACTION"}->[0],"LOAD");
441                            for (my $j=0; $j < @{$row->{MEDIA}};$j++) {
442                                    my $class = "Active <=>";
443                                    if ($row->{CLASS}->[$j] eq "Positive") {
444                                            $class = "Essential =>";
445                                    } elsif ($row->{CLASS}->[$j] eq "Negative") {
446                                            $class = "Essential <=";
447                                    } elsif ($row->{CLASS}->[$j] eq "Blocked") {
448                                            $class = "Inactive";
449                                    } elsif ($row->{CLASS}->[$j] eq "Positive variable") {
450                                            $class = "Active =>";
451                                    } elsif ($row->{CLASS}->[$j] eq "Negative variable") {
452                                            $class = "Active <=";
453                                    } elsif ($row->{CLASS}->[$j] eq "Variable") {
454                                            $class = "Active <=>";
455                                    } elsif ($row->{CLASS}->[$j] eq "Dead") {
456                                            $class = "Dead";
457                                    }
458                                    push(@{$rxnRow->{PREDICTIONS}},$row->{MEDIA}->[$j].":".$class);
459                            }
460                    }
461          }          }
462    
463          return $self->{_reaction_data};          return $self->{_reaction_data};
464  }  }
465    
466    =head3 feature_table
467    Definition:
468            FIGMODELTable = FIGMODELmodel->feature_table();
469    Description:
470            Returns FIGMODELTable with the feature list for the model
471    =cut
472    sub feature_table {
473            my ($self) = @_;
474    
475            if (!defined($self->{_feature_data})) {
476                    #Getting the genome feature list
477                    my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome());
478                    if (!defined($FeatureTable)) {
479                            print STDERR "FIGMODELmodel:feature_table:Could not get features for genome ".$self->genome()." in database!";
480                            return undef;
481                    }
482                    #Getting the reaction table for the model
483                    my $rxnTable = $self->reaction_table();
484                    if (!defined($rxnTable)) {
485                    print STDERR "FIGMODELmodel:feature_table:Could not get reaction table for model ".$self->id()." in database!";
486                            return undef;
487                }
488                    #Cloning the feature table
489                    $self->{_feature_data} = $FeatureTable->clone_table_def();
490                    $self->{_feature_data}->add_headings(($self->id()."REACTIONS",$self->id()."PREDICTIONS"));
491                for (my $i=0; $i < $rxnTable->size(); $i++) {
492                    my $Row = $rxnTable->get_row($i);
493                    if (defined($Row) && defined($Row->{"ASSOCIATED PEG"})) {
494                            foreach my $GeneSet (@{$Row->{"ASSOCIATED PEG"}}) {
495                                    my $temp = $GeneSet;
496                                    $temp =~ s/\+/|/g;
497                                    $temp =~ s/\sAND\s/|/gi;
498                                    $temp =~ s/\sOR\s/|/gi;
499                                    $temp =~ s/[\(\)\s]//g;
500                                    my @GeneList = split(/\|/,$temp);
501                                    foreach my $Gene (@GeneList) {
502                                    my $FeatureRow = $self->{_feature_data}->get_row_by_key("fig|".$self->genome().".".$Gene,"ID");
503                                                    if (!defined($FeatureRow)) {
504                                                            $FeatureRow = $FeatureTable->get_row_by_key("fig|".$self->genome().".".$Gene,"ID");
505                                                            if (defined($FeatureRow)) {
506                                                                    $self->{_feature_data}->add_row($FeatureRow);
507                                                            }
508                                                    }
509                                                    if (defined($FeatureRow)) {
510                                                            $self->{_feature_data}->add_data($FeatureRow,$self->id()."REACTIONS",$Row->{"LOAD"}->[0],1);
511                                                    }
512                                    }
513                            }
514                    }
515                }
516                #Loading predictions
517                my @files = glob($self->directory()."EssentialGenes-".$self->id()."-*");
518                foreach my $file (@files) {
519                    if ($file =~ m/\-([^\-]+)\.tbl/) {
520                            my $media = $1;
521                            my $list = $self->figmodel()->database()->load_single_column_file($file,"");
522                            my $hash;
523                            foreach my $gene (@{$list}) {
524                                    $hash->{$gene} = 1;
525                            }
526                            for (my $i=0; $i < $self->{_feature_data}->size(); $i++) {
527                                    my $Row = $self->{_feature_data}->get_row($i);
528                                    if ($Row->{ID}->[0] =~ m/(peg\.\d+)/) {
529                                            my $gene = $1;
530                                            if (defined($hash->{$gene})) {
531                                                    push(@{$Row->{$self->id()."PREDICTIONS"}},$media.":essential");
532                                            } else {
533                                                    push(@{$Row->{$self->id()."PREDICTIONS"}},$media.":nonessential");
534                                            }
535                                    }
536                            }
537                    }
538                }
539            }
540            return $self->{_feature_data};
541    }
542    
543  =head3 reaction_class_table  =head3 reaction_class_table
544  Definition:  Definition:
545          FIGMODELTable = FIGMODELmodel->reaction_class_table();          FIGMODELTable = FIGMODELmodel->reaction_class_table();
# Line 438  Line 580 
580          return $self->{_compound_class_table};          return $self->{_compound_class_table};
581  }  }
582    
583    =head3 get_essential_genes
584    Definition:
585            [string::peg ID] = FIGMODELmodel->get_essential_genes(string::media condition);
586    Description:
587            Returns an reference to an array of the predicted essential genes during growth in the input media condition
588    =cut
589    sub get_essential_genes {
590            my ($self,$media) = @_;
591    
592            if (!defined($media)) {
593                    $media = "Complete";
594            }
595            if (!defined($self->{_essential_genes}->{$media})) {
596                    $self->{_essential_genes}->{$media} = undef;
597                    if (-e $self->directory()."EssentialGenes-".$self->id().$self->selected_version()."-".$media.".tbl") {
598                            $self->{_essential_genes}->{$media} = $self->figmodel()->database()->load_single_column_file($self->directory()."EssentialGenes-".$self->id().$self->selected_version()."-".$media.".tbl","");
599                    }
600            }
601    
602            return $self->{_essential_genes}->{$media};
603    }
604    
605  =head3 compound_table  =head3 compound_table
606  Definition:  Definition:
607          FIGMODELTable = FIGMODELmodel->compound_table();          FIGMODELTable = FIGMODELmodel->compound_table();
# Line 456  Line 620 
620    
621  =head3 get_compound_data  =head3 get_compound_data
622  Definition:  Definition:
623          string = FIGMODELmodel->get_compound_data(string::compound ID);          {string:key=>[string]:values} = FIGMODELmodel->get_compound_data(string::compound ID);
624  Description:  Description:
625          Returns model compound data          Returns model compound data
626  =cut  =cut
# Line 467  Line 631 
631                  return undef;                  return undef;
632          }          }
633          if ($compound =~ m/^\d+$/) {          if ($compound =~ m/^\d+$/) {
634                  $self->compound_table()->get_row($compound);                  return $self->compound_table()->get_row($compound);
635          }          }
636          return $self->compound_table()->get_row_by_key($compound,"DATABASE");          return $self->compound_table()->get_row_by_key($compound,"DATABASE");
637  }  }
638    
639    =head3 get_feature_data
640    Definition:
641            {string:key=>[string]:values} = FIGMODELmodel->get_feature_data(string::feature ID);
642    Description:
643            Returns model feature data
644    =cut
645    sub get_feature_data {
646            my ($self,$feature) = @_;
647            if (!defined($self->feature_table())) {
648                    return undef;
649            }
650            if ($feature =~ m/^\d+$/) {
651                    return $self->feature_table()->get_row($feature);
652            }
653            if ($feature =~ m/(peg\.\d+)/) {
654                    $feature = $1;
655            }
656            return $self->feature_table()->get_row_by_key("fig|".$self->genome().".".$feature,"ID");
657    }
658    
659  =head3 genome  =head3 genome
660  Definition:  Definition:
661          string = FIGMODELmodel->genome();          string = FIGMODELmodel->genome();
# Line 480  Line 664 
664  =cut  =cut
665  sub genome {  sub genome {
666          my ($self) = @_;          my ($self) = @_;
667          return $self->{_data}->{genome}->[0];          return $self->{_data}->genome();
668  }  }
669    
670  =head3 source  =head3 source
# Line 491  Line 675 
675  =cut  =cut
676  sub source {  sub source {
677          my ($self) = @_;          my ($self) = @_;
678          return $self->{_data}->{source}->[0];          return $self->{_data}->source();
679  }  }
680    
681  =head3 rights  =head3 rights
# Line 502  Line 686 
686  =cut  =cut
687  sub rights {  sub rights {
688          my ($self,$username) = @_;          my ($self,$username) = @_;
   
689          if ($self->public()) {          if ($self->public()) {
690                  return 1;                  return 1;
691          }          }
# Line 510  Line 693 
693                  return 0;                  return 0;
694          }          }
695          if (!defined($self->{_userrights}->{$username})) {          if (!defined($self->{_userrights}->{$username})) {
696                  if (defined($self->{_data}->{master})) {                  $self->{_userrights}->{$self->{_data}->owner()} = 1;
697                          for (my $i=0; $i < @{$self->{_data}->{master}};$i++) {                  my @users = split(/\|/,$self->{_data}->users());
698                                  $self->{_userrights}->{$self->{_data}->{master}->[$i]} = 1;                  for (my $i=0; $i < @users;$i++) {
699                          }                          $self->{_userrights}->{$users[$i]} = 1;
                 }  
                 if (defined($self->{_data}->{users})) {  
                         for (my $i=0; $i < @{$self->{_data}->{users}};$i++) {  
                                 $self->{_userrights}->{$self->{_data}->{users}->[$i]} = 1;  
                         }  
700                  }                  }
701          }          }
702          return $self->{_userrights}->{$username};          return $self->{_userrights}->{$username};
# Line 532  Line 710 
710  =cut  =cut
711  sub public {  sub public {
712          my ($self) = @_;          my ($self) = @_;
713            if ($self->{_data}->users() eq "all") {
714          if (!defined($self->{_public})) {                  return 1;
                 $self->{_public} = 0;  
                 if (defined($self->{_data}->{users}->[0]) && $self->{_data}->{users}->[0] eq "all") {  
                         $self->{_public} = 1;  
                 }  
715          }          }
716          return $self->{_public};          return 0;
717  }  }
718    
719  =head3 directory  =head3 directory
# Line 587  Line 761 
761          return $self->directory().$self->id().$self->selected_version().".txt";          return $self->directory().$self->id().$self->selected_version().".txt";
762  }  }
763    
 =head3 set_metagenome_stats  
 Definition:  
         string = FIGMODELmodel->set_metagenome_stats();  
 Description:  
         Sets the values of many model stats for a metagenome  
 =cut  
 sub set_metagenome_stats {  
         my ($self) = @_;  
   
         $self->{_total_compounds} = 0;  
         if (defined($self->compound_table())) {  
                 $self->{_total_compounds} = $self->compound_table()->size();  
         }  
         $self->{_gene_reactions} = 0;  
         $self->{_gapfilling_reactions} = 0;  
         $self->{_model_genes} = 0;  
         $self->{_total_reactions} = 0;  
         if (defined($self->reaction_table())) {  
                 $self->{_total_reactions} = $self->reaction_table()->size();  
                 my $tbl = $self->reaction_table();  
                 my $spontaneous = 0;  
                 for (my $i=0; $i < $tbl->size(); $i++) {  
                         my $row = $tbl->get_row($i);  
                         if (!defined($row->{"ASSOCIATED PEG"}->[0]) || $row->{"ASSOCIATED PEG"}->[0] !~ m/peg/) {  
                                 if ($row->{"ASSOCIATED PEG"}->[0] =~ m/SPONTANEOUS/) {  
                                         $spontaneous++;  
                                 } else {  
                                         $self->{_gapfilling_reactions}++;  
                                 }  
                         } else {  
                                 for (my $j=0; $j < @{$row->{"CONFIDENCE"}}; $j++) {  
                                         my @ecores = split(/;/,$row->{"CONFIDENCE"}->[$j]);  
                                         $self->{_model_genes} += @ecores;  
                                 }  
                         }  
                 }  
                 $self->{_gene_reactions} = $tbl->size() - $spontaneous - $self->{_gapfilling_reactions};  
         }  
 }  
   
764  =head3 version  =head3 version
765  Definition:  Definition:
766          string = FIGMODELmodel->version();          string = FIGMODELmodel->version();
# Line 638  Line 772 
772    
773          if (!defined($self->{_version})) {          if (!defined($self->{_version})) {
774                  if (!defined($self->{_selectedversion})) {                  if (!defined($self->{_selectedversion})) {
775                          if (defined($self->stats())) {                          $self->{_version} = "V".$self->{_data}->version().".".$self->{_data}->autocompleteVersion();
                                 $self->{_version} = "V".$self->stats()->{"Version"}->[0].".".$self->stats()->{"Gap fill version"}->[0];  
                         }  
776                  } else {                  } else {
777                          $self->{_version} = $self->{_selectedversion};                          $self->{_version} = $self->{_selectedversion};
778                  }                  }
# Line 671  Line 803 
803  =cut  =cut
804  sub modification_time {  sub modification_time {
805          my ($self) = @_;          my ($self) = @_;
806          if (!defined($self->{_modification_time})) {          return $self->{_data}->modificationDate();
                 my $stats = $self->stats();  
                 if (defined($stats)) {  
                         $self->{_modification_time} = 0;  
                         if (defined($stats->{"Build date"}->[0]) && $self->{_modification_time} < $stats->{"Build date"}->[0]) {  
                                 $self->{_modification_time} = $stats->{"Build date"}->[0];  
                         } elsif (defined($stats->{"Gap fill date"}->[0]) && $self->{_modification_time} < $stats->{"Gap fill date"}->[0]) {  
                                 $self->{_modification_time} = $stats->{"Gap fill date"}->[0];  
                         }  
                 } else {  
                         $self->{_modification_time} = 0;  
                 }  
         }  
         return $self->{_modification_time};  
807  }  }
808    
809  =head3 gene_reactions  =head3 gene_reactions
# Line 695  Line 814 
814  =cut  =cut
815  sub gene_reactions {  sub gene_reactions {
816          my ($self) = @_;          my ($self) = @_;
817            return ($self->{_data}->reactions() - $self->{_data}->autoCompleteReactions() - $self->{_data}->spontaneousReactions() - $self->{_data}->gapFillReactions());
         if (!defined($self->{_gene_reactions})) {  
                 if ($self->source() =~ /^MGRAST/) {  
                         $self->set_metagenome_stats();  
                 } elsif (defined($self->stats())) {  
                         $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];  
                 }  
         }  
         return $self->{_gene_reactions};  
818  }  }
819    
820  =head3 total_compounds  =head3 total_compounds
# Line 714  Line 825 
825  =cut  =cut
826  sub total_compounds {  sub total_compounds {
827          my ($self) = @_;          my ($self) = @_;
828            return $self->{_data}->compounds();
         if (!defined($self->{_total_compounds})) {  
                 if ($self->source() =~ /^MGRAST/) {  
                         $self->set_metagenome_stats();  
                 } elsif (defined($self->stats())) {  
                         $self->{_total_compounds} = $self->stats()->{'Metabolites'}->[0];  
                 }  
         }  
         return $self->{_total_compounds};  
829  }  }
830    
831  =head3 gapfilling_reactions  =head3 gapfilling_reactions
# Line 733  Line 836 
836  =cut  =cut
837  sub gapfilling_reactions {  sub gapfilling_reactions {
838          my ($self) = @_;          my ($self) = @_;
839            return ($self->{_data}->autoCompleteReactions()+$self->{_data}->gapFillReactions());
         if (!defined($self->{_gapfilling_reactions})) {  
                 if ($self->source() =~ /^MGRAST/) {  
                         $self->set_metagenome_stats();  
                 } elsif (defined($self->stats())) {  
                         $self->{_gapfilling_reactions} = $self->stats()->{'Gap filling reactions'}->[0];  
                 }  
         }  
         return $self->{_gapfilling_reactions};  
840  }  }
841    
842  =head3 total_reactions  =head3 total_reactions
# Line 752  Line 847 
847  =cut  =cut
848  sub total_reactions {  sub total_reactions {
849          my ($self) = @_;          my ($self) = @_;
850            return $self->{_data}->reactions();
         if (!defined($self->{_total_reactions})) {  
                 if ($self->source() =~ /^MGRAST/) {  
                         $self->set_metagenome_stats();  
                 } elsif (defined($self->stats())) {  
                         $self->{_total_reactions} = $self->stats()->{'Number of reactions'}->[0];  
                 }  
         }  
         return $self->{_total_reactions};  
851  }  }
852    
853  =head3 model_genes  =head3 model_genes
# Line 771  Line 858 
858  =cut  =cut
859  sub model_genes {  sub model_genes {
860          my ($self) = @_;          my ($self) = @_;
861            return $self->{_data}->associatedGenes();
         if (!defined($self->{_model_genes})) {  
                 if ($self->source() =~ /^MGRAST/) {  
                         $self->set_metagenome_stats();  
                 } elsif (defined($self->stats())) {  
                         $self->{_model_genes} = $self->stats()->{'Genes with reactions'}->[0];  
                 }  
         }  
         return $self->{_model_genes};  
862  }  }
863    
864  =head3 class  =head3 class
# Line 790  Line 869 
869  =cut  =cut
870  sub class {  sub class {
871          my ($self) = @_;          my ($self) = @_;
872            return $self->{_data}->cellwalltype();
         if (!defined($self->{_class})) {  
                 if ($self->source() =~ /^MGRAST/) {  
                         $self->{_class} = "Other";  
                 } elsif (defined($self->stats())) {  
                         $self->{_class} = $self->stats()->{Class}->[0];  
                 }  
         }  
         return $self->{_class};  
873  }  }
874    
875  =head3 taxonomy  =head3 taxonomy
# Line 866  Line 937 
937                          if (defined($self->mgdata())) {                          if (defined($self->mgdata())) {
938                                  $self->{_genome_genes} = $self->mgdata()->genome_contig_count;                                  $self->{_genome_genes} = $self->mgdata()->genome_contig_count;
939                          }                          }
940                  } elsif (defined($self->stats())) {                  } else {
941                          $self->{_genome_genes} = $self->stats()->{'Total genes'}->[0];                          $self->{_genome_genes} = $self->figmodel()->get_genome_stats($self->genome())->{"TOTAL GENES"}->[0];
942                  }                  }
943          }          }
944    
# Line 892  Line 963 
963          #Checking that the table is defined and the output file exists          #Checking that the table is defined and the output file exists
964          if (defined($result) && defined($result->get_row(0)->{"ESSENTIALGENES"})) {          if (defined($result) && defined($result->get_row(0)->{"ESSENTIALGENES"})) {
965                  $self->figmodel()->database()->print_array_to_file($self->directory()."EssentialGenes-".$self->id()."-".$Media.".tbl",[join("\n",@{$result->get_row(0)->{"ESSENTIALGENES"}})]);                  $self->figmodel()->database()->print_array_to_file($self->directory()."EssentialGenes-".$self->id()."-".$Media.".tbl",[join("\n",@{$result->get_row(0)->{"ESSENTIALGENES"}})]);
966            } else {
967                    $self->figmodel()->error_message("FIGMODELmodel:run_default_model_predictions:could not identify essential reactions for model ".$self->id().$self->selected_version().".");
968                    return $self->figmodel()->fail();
969          }          }
970    
971          #Classifying reactions and compounds          #Classifying reactions and compounds
972          my $tbl = $self->classify_model_reactions($Media);          my $tbl = $self->classify_model_reactions($Media);
973            if (!defined($tbl)) {
974                    $self->figmodel()->error_message("FIGMODELmodel:run_default_model_predictions:could not classify reactions for model ".$self->id().$self->selected_version().".");
975                    return $self->figmodel()->fail();
976            }
977          $tbl->save();          $tbl->save();
978    
979          return $self->figmodel()->success();          return $self->figmodel()->success();
980  }  }
981    
 =head3 save_obsolete_stats  
 Definition:  
         FIGMODELmodel->save_obsolete_stats();  
 Description:  
 =cut  
 sub save_obsolete_stats {  
         my ($self) = @_;  
   
         #checking if stats exists  
         my $stats = $self->stats();  
         if (defined($stats)) {  
                 $stats->{"Model ID"}->[0] = $self->id()."V".$stats->{"Version"}->[0].".".$stats->{"Gap fill version"}->[0];  
                 $self->figmodel()->database()->update_row("OBSOLETE MODEL STATS",$stats,"Model ID");  
                 $stats->{"Model ID"}->[0] = $self->id();  
         }  
 }  
   
982  =head3 update_stats_for_gap_filling  =head3 update_stats_for_gap_filling
983  Definition:  Definition:
984          {string => [string]} = FIGMODELmodel->update_stats_for_gap_filling(int::gapfill time);          {string => [string]} = FIGMODELmodel->update_stats_for_gap_filling(int::gapfill time);
# Line 925  Line 986 
986  =cut  =cut
987  sub update_stats_for_gap_filling {  sub update_stats_for_gap_filling {
988          my ($self,$gapfilltime) = @_;          my ($self,$gapfilltime) = @_;
989            $self->{_data}->autoCompleteTime($gapfilltime);
990          #preserving the stats for the now obselete model          $self->{_data}->autocompleteDate(time());
991          $self->save_obsolete_stats();          $self->{_data}->modificationDate(time());
992          my $stats = $self->update_model_stats(0);          my $version = $self->{_data}->autocompleteVersion();
993          $stats->{"Gap filling time"}->[0] = $gapfilltime;          $self->{_data}->autocompleteVersion($version+1);
         $stats->{"Gap fill date"}->[0] = time();  
         if (!defined($stats->{"Gap fill version"}->[0])) {  
                 $stats->{"Gap fill version"}->[0] = 0;  
         }  
         $stats->{"Gap fill version"}->[0]++;  
   
         #Updating the stats stored in the table  
         $self->figmodel()->database()->update_row("MODEL STATS",$stats,"Model ID");  
         return $stats;  
994  }  }
995    
996  =head3 update_stats_for_build  =head3 update_stats_for_build
# Line 948  Line 1000 
1000  =cut  =cut
1001  sub update_stats_for_build {  sub update_stats_for_build {
1002          my ($self) = @_;          my ($self) = @_;
1003            $self->{_data}->builtDate(time());
1004          #preserving the stats for the now obselete model          $self->{_data}->modificationDate(time());
1005          $self->save_obsolete_stats();          my $version = $self->{_data}->version();
1006          my $stats = $self->update_model_stats(0);          $self->{_data}->version($version+1);
         $stats->{"Build date"}->[0] = time();  
         if (!defined($stats->{"Version"}->[0])) {  
                 $stats->{"Version"}->[0] = 0;  
         }  
         $stats->{"Version"}->[0]++;  
   
         #Updating the stats stored in the table  
         $self->figmodel()->database()->update_row("MODEL STATS",$stats,"Model ID");  
         return $stats;  
1007  }  }
1008    
1009  =head3 update_model_stats  =head3 update_model_stats
# Line 979  Line 1022 
1022          }          }
1023          my $cpdtbl = $self->compound_table();          my $cpdtbl = $self->compound_table();
1024    
1025          #Creating empty status row          #Calculating all necessary stats
         my $CurrentStats = {"Genes with reactions" => [0],  
                                                  "Metabolites" => [$cpdtbl->size()],  
                                                  "Growmatch reactions" => [0],  
                                                  "Spontaneous" => [0],  
                                                  "Biolog gap filling reactions" => [0],  
                                                  "Number of reactions" => [$rxntbl->size()],  
                                                  "Transport reaction"=>[0],  
                                                  "Gap filling reactions" => [0],  
                                                  "Model ID" => [$self->id()],  
                                                  "Subsystem genes with reactions" => [0],  
                                                  "Nonsubsystem genes with reactions" => [0],  
                                                  Version => [0],  
                                                  "Gap fill version" => [0],  
                                                  Source => [$self->source()],  
                                                  "Genes with one reaction" => [0],  
                                                  "Genome ID" => [$self->genome()]};  
   
         my $genomestats = $self->figmodel()->get_genome_stats($self->genome());  
         if (defined($genomestats)) {  
                 $CurrentStats->{SOURCE}->[0] = $genomestats->{SOURCE}->[0];  
                 $CurrentStats->{"Organism name"}->[0] = $genomestats->{NAME}->[0];  
                 $CurrentStats->{"Total genes"}->[0] = $genomestats->{"TOTAL GENES"}->[0];  
                 $CurrentStats->{"Gram positive genes"}->[0] = $genomestats->{"GRAM POSITIVE GENES"}->[0];  
                 $CurrentStats->{"Gram negative genes"}->[0] = $genomestats->{"GRAM NEGATIVE GENES"}->[0];  
                 $CurrentStats->{"Class"}->[0] = $genomestats->{CLASS}->[0];  
                 $CurrentStats->{"Genes with functions"}->[0] = $genomestats->{"GENES WITH FUNCTIONS"}->[0];  
                 $CurrentStats->{"Subsystem genes"}->[0] = $genomestats->{"SUBSYSTEM GENES"}->[0];  
                 $CurrentStats->{"Nonsubsystem genes"}->[0] = $genomestats->{"NON SUBSYSTEM GENES"}->[0];  
                 if (defined($genomestats->{TAXONOMY})) {  
                         $CurrentStats->{"Taxonomy 0"}->[0] = $genomestats->{TAXONOMY}->[0];  
                         $CurrentStats->{"Taxonomy 1"}->[0] = $genomestats->{TAXONOMY}->[1];  
                         $CurrentStats->{"Taxonomy 2"}->[0] = $genomestats->{TAXONOMY}->[2];  
                         $CurrentStats->{"Taxonomy 3"}->[0] = $genomestats->{TAXONOMY}->[3];  
                         $CurrentStats->{"Taxonomy 4"}->[0] = $genomestats->{TAXONOMY}->[4];  
                         $CurrentStats->{"Taxonomy 5"}->[0] = $genomestats->{TAXONOMY}->[5];  
                 }  
         }  
   
         #Transfering build, version, and gap fill data from existing stats  
         if (defined($self->stats())) {  
                 $CurrentStats->{Version}->[0] = $self->stats()->{Version}->[0];  
                 $CurrentStats->{"Gap fill version"}->[0] = $self->stats()->{"Gap fill version"}->[0];  
                 $CurrentStats->{"Gap filling time"}->[0] = $self->stats()->{"Gap filling time"}->[0];  
                 $CurrentStats->{"Gap fill date"}->[0] = $self->stats()->{"Gap fill date"}->[0];  
                 $CurrentStats->{"Build date"}->[0] = $self->stats()->{"Build date"}->[0];  
         }  
   
1026          my %GeneHash;          my %GeneHash;
1027          my %NonpegHash;          my %NonpegHash;
1028          my %CompoundHash;          my %CompoundHash;
1029            my $spontaneousReactions = 0;
1030            my $gapFillReactions = 0;
1031            my $biologReactions = 0;
1032            my $transporters = 0;
1033            my $autoCompleteReactions = 0;
1034            my $associatedSubsystemGenes = 0;
1035          for (my $i=0; $i < $rxntbl->size(); $i++) {          for (my $i=0; $i < $rxntbl->size(); $i++) {
1036                  my $Row = $rxntbl->get_row($i);                  my $Row = $rxntbl->get_row($i);
1037                  if (defined($Row) && defined($Row->{"ASSOCIATED PEG"})) {                  if (defined($Row) && defined($Row->{"ASSOCIATED PEG"})) {
# Line 1037  Line 1039 
1039                          if (defined($ReactionRow->{"EQUATION"}->[0])) {                          if (defined($ReactionRow->{"EQUATION"}->[0])) {
1040                                  #Checking for extracellular metabolites which indicate that this is a transporter                                  #Checking for extracellular metabolites which indicate that this is a transporter
1041                                  if ($ReactionRow->{"EQUATION"}->[0] =~ m/\[e\]/) {                                  if ($ReactionRow->{"EQUATION"}->[0] =~ m/\[e\]/) {
1042                                          $CurrentStats->{"Transport reaction"}->[0]++;                                          $transporters++;
1043                                  }                                  }
1044                          }                          }
1045                          #Identifying spontaneous/biolog/gapfilling/gene associated reactions                          #Identifying spontaneous/biolog/gapfilling/gene associated reactions
1046                          if ($Row->{"ASSOCIATED PEG"}->[0] =~ m/BIOLOG/i) {                          if ($Row->{"ASSOCIATED PEG"}->[0] =~ m/BIOLOG/i) {
1047                                  $CurrentStats->{"Biolog gap filling reactions"}->[0]++;                                  $biologReactions++;
1048                            } elsif ($Row->{"ASSOCIATED PEG"}->[0] =~ m/GROW/i) {
1049                                    $gapFillReactions++;
1050                          } elsif ($Row->{"ASSOCIATED PEG"}->[0] =~ m/SPONTANEOUS/i) {                          } elsif ($Row->{"ASSOCIATED PEG"}->[0] =~ m/SPONTANEOUS/i) {
1051                                  $CurrentStats->{"Spontaneous"}->[0]++;                                  $spontaneousReactions++;
1052                          } elsif ($Row->{"ASSOCIATED PEG"}->[0] =~ m/GAP/ || $Row->{"ASSOCIATED PEG"}->[0] =~ m/UNIVERSAL/i || $Row->{"ASSOCIATED PEG"}->[0] =~ m/UNKNOWN/i) {                          } elsif ($Row->{"ASSOCIATED PEG"}->[0] =~ m/GAP/ || $Row->{"ASSOCIATED PEG"}->[0] =~ m/UNIVERSAL/i || $Row->{"ASSOCIATED PEG"}->[0] =~ m/UNKNOWN/i) {
1053                                  $CurrentStats->{"Gap filling reactions"}->[0]++;                                  $autoCompleteReactions++;
1054                          } else {                          } else {
1055                                  foreach my $GeneSet (@{$Row->{"ASSOCIATED PEG"}}) {                                  foreach my $GeneSet (@{$Row->{"ASSOCIATED PEG"}}) {
1056                                          my @GeneList = split(/\+/,$GeneSet);                                          $_ = $GeneSet;
1057                                            my @GeneList = /(peg\.\d+)/g;
1058                                          foreach my $Gene (@GeneList) {                                          foreach my $Gene (@GeneList) {
1059                                                  if ($Gene =~ m/(peg\.\d+)/) {                                                  if ($Gene =~ m/(peg\.\d+)/) {
1060                                                          $GeneHash{$1} = 1;                                                          $GeneHash{$1} = 1;
# Line 1063  Line 1068 
1068          }          }
1069          my @genes = keys(%GeneHash);          my @genes = keys(%GeneHash);
1070          my @othergenes = keys(%NonpegHash);          my @othergenes = keys(%NonpegHash);
         $CurrentStats->{"Genes with reactions"}->[0] = @genes + @othergenes;  
1071    
1072          #Updating the stats stored in the table          #Setting the reaction count
1073          $self->figmodel()->database()->update_row("MODEL STATS",$CurrentStats,"Model ID");          $self->{_data}->reactions($rxntbl->size());
1074          $self->{_stats} = $CurrentStats;          #Setting the metabolite count
1075          return $CurrentStats;          $self->{_data}->compounds($rxntbl->size());
1076            #Setting the gene count
1077            my $geneCount = @genes + @othergenes;
1078            $self->{_data}->associatedGenes($geneCount);
1079            #Setting remaining stats
1080            $self->{_data}->spontaneousReactions($spontaneousReactions);
1081            $self->{_data}->gapFillReactions($gapFillReactions);
1082            $self->{_data}->biologReactions($biologReactions);
1083            $self->{_data}->transporters($transporters);
1084            $self->{_data}->autoCompleteReactions($autoCompleteReactions);
1085            $self->{_data}->associatedSubsystemGenes($associatedSubsystemGenes);
1086            #Setting the model class
1087            my $class = "";
1088            for (my $i=0; $i < @{$self->figmodel()->config("class list")}; $i++) {
1089                    if (defined($self->figmodel()->config($self->figmodel()->config("class list")->[$i]))) {
1090                            if (defined($self->figmodel()->config($self->figmodel()->config("class list")->[$i])->{$self->id()})) {
1091                                    $class = $self->figmodel()->config("class list")->[$i];
1092                                    last;
1093                            }
1094                            if ($class eq "" && defined($self->figmodel()->config($self->figmodel()->config("class list")->[$i])->{$self->genome()})) {
1095                                    $class = $self->figmodel()->config("class list")->[$i];
1096                            }
1097                    }
1098            }
1099            if ($class eq "") {
1100                    $class = $self->figmodel()->get_genome_stats($self->genome())->{CLASS}->[0];
1101            }
1102            if ($class eq "") {
1103                    $class = "unknown";
1104            }
1105            $self->{_data}->cellwalltype($class);
1106  }  }
1107    
1108  =head3 GapFillModel  =head3 GapFillModel
# Line 1119  Line 1153 
1153                          }                          }
1154                  }                  }
1155                  $MediaTable->save($self->config("Media directory")->[0].$UniqueFilename."TestMedia.txt");                  $MediaTable->save($self->config("Media directory")->[0].$UniqueFilename."TestMedia.txt");
1156                  system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$UniqueFilename."TestMedia",["GapFilling"],{"Default max drain flux" => 0},"GapFill".$self->id().".log",undef));                  system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$UniqueFilename."TestMedia",["GapFilling"],{"Default max drain flux" => 0,"Reactions to knockout" => $self->config("permanently knocked out reactions")->[0]},"GapFill".$self->id().".log",undef));
1157                  unlink($self->config("Media directory")->[0].$UniqueFilename."TestMedia.txt");                  unlink($self->config("Media directory")->[0].$UniqueFilename."TestMedia.txt");
1158          } else {          } else {
1159                  system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),undef,["GapFilling"],undef,"GapFill".$self->id().".log",undef));                  system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),undef,["GapFilling"],{"Reactions to knockout" => $self->config("permanently knocked out reactions")->[0]},"GapFill".$self->id().".log",undef));
1160          }          }
1161    
1162          #Parse the solutions file for the model and get the reaction list from it          #Parse the solutions file for the model and get the reaction list from it
# Line 1140  Line 1174 
1174          for (my $i=($SolutionData->size()-1); $i >=0; $i--) {          for (my $i=($SolutionData->size()-1); $i >=0; $i--) {
1175                  if (defined($SolutionData->get_row($i)->{"Notes"}) && $SolutionData->get_row($i)->{"Notes"}->[0] =~ m/^Recursive/) {                  if (defined($SolutionData->get_row($i)->{"Notes"}) && $SolutionData->get_row($i)->{"Notes"}->[0] =~ m/^Recursive/) {
1176                          my $AllSolutions = substr($SolutionData->get_row($i)->{"Notes"}->[0],15);                          my $AllSolutions = substr($SolutionData->get_row($i)->{"Notes"}->[0],15);
                         print "Solution:".$AllSolutions."\n";  
1177                          my @TempThree = split(/\|/,$AllSolutions);                          my @TempThree = split(/\|/,$AllSolutions);
1178                          if (@TempThree > 0 && $TempThree[0] =~ m/.+:(.+)/) {                          if (@TempThree > 0 && $TempThree[0] =~ m/.+:(.+)/) {
1179                                  my @TempFour = split(/,/,$1);                                  my @TempFour = split(/,/,$1);
# Line 1163  Line 1196 
1196                                                  push(@{$ReactionList},$ID);                                                  push(@{$ReactionList},$ID);
1197                                          }                                          }
1198                                  }                                  }
                                 print "Solution:".$TempThree[0]."\n";  
1199                                  $self->figmodel()->IntegrateGrowMatchSolution($self->id(),undef,$ReactionList,$DirectionList,"GAP FILLING",0,1);                                  $self->figmodel()->IntegrateGrowMatchSolution($self->id(),undef,$ReactionList,$DirectionList,"GAP FILLING",0,1);
1200                          }                          }
1201                          last;                          last;
# Line 1171  Line 1203 
1203          }          }
1204          #Updating model stats with gap filling results          #Updating model stats with gap filling results
1205          my $ElapsedTime = time() - $StartTime;          my $ElapsedTime = time() - $StartTime;
1206          $self->figmodel()->ClearDBModel($self->id(),1);          $self->figmodel()->database()->ClearDBModel($self->id(),1);
1207          #Determining why each gap filling reaction was added          #Determining why each gap filling reaction was added
1208          $self->figmodel()->IdentifyDependancyOfGapFillingReactions($self->id(),$Media);          $self->figmodel()->IdentifyDependancyOfGapFillingReactions($self->id(),$Media);
1209          if (!defined($donotclear) || $donotclear != 1) {          if (!defined($donotclear) || $donotclear != 1) {
# Line 1191  Line 1223 
1223          return $self->success();          return $self->success();
1224  }  }
1225    
1226    =head3 GapGenModel
1227    Definition:
1228            FIGMODELmodel->GapGenModel();
1229    Description:
1230            Runs the gap generation algorithm to correct a single false positive prediction. Results are loaded into a table.
1231    =cut
1232    
1233    sub GapGenModel {
1234            my ($self,$Media,$KOList,$NoKOList,$Experiment,$SolutionLimit) = @_;
1235    
1236            #Enforcing nonoptional arguments
1237            if (!defined($Media)) {
1238                    return undef;
1239            }
1240            if (!defined($KOList)) {
1241                    $KOList->[0] = "none";
1242            }
1243            if (!defined($NoKOList)) {
1244                    $NoKOList->[0] = "none";
1245            }
1246            if (!defined($Experiment)) {
1247                    $Experiment= "ReactionKO";
1248            }
1249            if (!defined($SolutionLimit)) {
1250                    $SolutionLimit = "10";
1251            }
1252    
1253            #Translating the KO lists into arrays
1254            if (ref($KOList) ne "ARRAY") {
1255                    my $temp = $KOList;
1256                    $KOList = ();
1257                    push(@{$KOList},split(/[,;]/,$temp));
1258            }
1259            my $noKOHash;
1260            if (defined($NoKOList) && ref($NoKOList) ne "ARRAY") {
1261                    my $temp = $NoKOList;
1262                    $NoKOList = ();
1263                    push(@{$NoKOList},split(/[,;]/,$temp));
1264                    foreach my $rxn (@{$NoKOList}) {
1265                            $noKOHash->{$rxn} = 1;
1266                    }
1267            }
1268    
1269            #Checking if solutions exist for the input parameters
1270            $self->aquireModelLock();
1271            my $tbl = $self->load_model_table("GapGenSolutions");
1272            my $solutionRow = $tbl->get_table_by_key($Experiment,"Experiment")->get_table_by_key($Media,"Media")->get_row_by_key(join(",",@{$KOList}),"KOlist");
1273            my $solutions;
1274            if (defined($solutionRow)) {
1275                    #Checking if any solutions conform to the no KO list
1276                    foreach my $solution (@{$solutionRow->{Solutions}}) {
1277                            my @reactions = split(/,/,$solution);
1278                            my $include = 1;
1279                            foreach my $rxn (@reactions) {
1280                                    if ($rxn =~ m/(rxn\d\d\d\d\d)/) {
1281                                            if (defined($noKOHash->{$1})) {
1282                                                    $include = 0;
1283                                            }
1284                                    }
1285                            }
1286                            if ($include == 1) {
1287                                    push(@{$solutions},$solution);
1288                            }
1289                    }
1290            } else {
1291                    $solutionRow = {Media => [$Media],Experiment => [$Experiment],KOlist => [join(",",@{$KOList})]};
1292                    $tbl->add_row($solutionRow);
1293                    $self->figmodel()->database()->save_table($tbl);
1294            }
1295            $self->releaseModelLock();
1296    
1297            #Returning solution list of solutions were found
1298            if (defined($solutions) && @{$solutions} > 0) {
1299                    return $solutions;
1300            }
1301    
1302            #Getting unique filename
1303            my $Filename = $self->figmodel()->filename();
1304    
1305            #Running the gap generation
1306            system($self->figmodel()->GenerateMFAToolkitCommandLineCall($Filename,$self->id().$self->selected_version(),$Media,["GapGeneration"],{"Recursive MILP solution limit" => $SolutionLimit ,"Reactions that should always be active" => join(";",@{$NoKOList}),"Reactions to knockout" => join(";",@{$KOList}),"Reactions that are always blocked" => "none"},"Gapgeneration-".$self->id().$self->selected_version()."-".$Filename.".log",undef,undef));
1307            my $ProblemReport = $self->figmodel()->LoadProblemReport($Filename);
1308            if (!defined($ProblemReport)) {
1309                    $self->figmodel()->error_message("FIGMODEL:GapGenerationAlgorithm;No problem report;".$Filename.";".$self->id().$self->selected_version().";".$Media.";".$KOList.";".$NoKOList);
1310                    return undef;
1311            }
1312    
1313            #Clearing the output folder and log file
1314            $self->figmodel()->clearing_output($Filename,"Gapgeneration-".$self->id().$self->selected_version()."-".$Filename.".log");
1315    
1316            #Saving the solution
1317            $self->aquireModelLock();
1318            $tbl = $self->load_model_table("GapGenSolutions");
1319            $solutionRow = $tbl->get_table_by_key($Experiment,"Experiment")->get_table_by_key($Media,"Media")->get_row_by_key(join(",",@{$KOList}),"KOlist");
1320            for (my $j=0; $j < $ProblemReport->size(); $j++) {
1321                    if ($ProblemReport->get_row($j)->{"Notes"}->[0] =~ m/^Recursive\sMILP\s([^)]+)/) {
1322                            my @SolutionList = split(/\|/,$1);
1323                            for (my $k=0; $k < @SolutionList; $k++) {
1324                                    if ($SolutionList[$k] =~ m/(\d+):(.+)/) {
1325                                            push(@{$solutionRow->{Solutions}},$2);
1326                                            push(@{$solutions},$2);
1327                                    }
1328                            }
1329                    }
1330            }
1331            $self->figmodel()->database()->save_table($tbl);
1332            $self->releaseModelLock();
1333    
1334            return $solutions;
1335    }
1336    
1337  =head3 datagapfill  =head3 datagapfill
1338  Definition:  Definition:
1339          success()/fail() = FIGMODELmodel->datagapfill();          success()/fail() = FIGMODELmodel->datagapfill();
# Line 1201  Line 1344 
1344          my ($self,$GapFillingRunSpecs,$TansferFileSuffix) = @_;          my ($self,$GapFillingRunSpecs,$TansferFileSuffix) = @_;
1345          my $UniqueFilename = $self->figmodel()->filename();          my $UniqueFilename = $self->figmodel()->filename();
1346          if (defined($GapFillingRunSpecs) && @{$GapFillingRunSpecs} > 0) {          if (defined($GapFillingRunSpecs) && @{$GapFillingRunSpecs} > 0) {
                 print "Gap filling specs:\n".join("\n",@{$GapFillingRunSpecs})."\n\n";  
1347                  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));                  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));
   
1348                  #Checking that the solution exists                  #Checking that the solution exists
1349                  if (!-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt") {                  if (!-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt") {
1350                          $self->figmodel()->error_message("FIGMODEL:GapFillingAlgorithm: Could not find MFA output file!");                          $self->figmodel()->error_message("FIGMODEL:GapFillingAlgorithm: Could not find MFA output file!");
# Line 1213  Line 1354 
1354                  my $GapFillResultTable = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt",";","",0,undef);                  my $GapFillResultTable = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt",";","",0,undef);
1355                  if (defined($TansferFileSuffix)) {                  if (defined($TansferFileSuffix)) {
1356                          system("cp ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt ".$self->directory().$self->id().$self->selected_version()."-".$TansferFileSuffix.".txt");                          system("cp ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt ".$self->directory().$self->id().$self->selected_version()."-".$TansferFileSuffix.".txt");
                         system("cp ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingReport.txt ".$self->directory().$self->id().$self->selected_version()."-".$TansferFileSuffix.".txt");  
1357                  }                  }
   
1358                  #If the system is not configured to preserve all logfiles, then the mfatoolkit output folder is deleted                  #If the system is not configured to preserve all logfiles, then the mfatoolkit output folder is deleted
1359                  $self->figmodel()->clearing_output($UniqueFilename,"GapFilling-".$self->id().$self->selected_version()."-".$UniqueFilename.".log");                  $self->figmodel()->clearing_output($UniqueFilename,"GapFilling-".$self->id().$self->selected_version()."-".$UniqueFilename.".log");
1360                  return $GapFillingRunSpecs;                  return $GapFillResultTable;
1361          }          }
1362          if (defined($TansferFileSuffix)) {          if (defined($TansferFileSuffix)) {
1363                  $self->figmodel()->database()->print_array_to_file($self->directory().$self->id().$self->selected_version()."-".$TansferFileSuffix.".txt",["Experiment;Solution index;Solution cost;Solution reactions"]);                  $self->figmodel()->database()->print_array_to_file($self->directory().$self->id().$self->selected_version()."-".$TansferFileSuffix.".txt",["Experiment;Solution index;Solution cost;Solution reactions"]);
# Line 1284  Line 1423 
1423                          push(@{$DirectionArray},$SolutionHash{$ReactionList[$k]});                          push(@{$DirectionArray},$SolutionHash{$ReactionList[$k]});
1424                  }                  }
1425                  print "Integrating solution!\n";                  print "Integrating solution!\n";
1426                  $self->IntegrateGrowMatchSolution($self->id().$self->selected_version(),$self->directory().$self->id().$TempVersion.".txt",$ReactionArray,$DirectionArray,"Gapfilling ".$GapFillResultTable->get_row($i)->{"Experiment"}->[0],1,1);                  $self->figmodel()->IntegrateGrowMatchSolution($self->id().$self->selected_version(),$self->directory().$self->id().$TempVersion.".txt",$ReactionArray,$DirectionArray,"Gapfilling ".$GapFillResultTable->get_row($i)->{"Experiment"}->[0],1,1);
1427                  my $model = $self->get_model($self->id().$TempVersion);                  $self->PrintModelLPFile();
                 $model->PrintModelLPFile();  
1428                  #Running the model against all available experimental data                  #Running the model against all available experimental data
1429                  print "Running test model!\n";                  print "Running test model!\n";
1430                  my ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector) = $self->RunAllStudiesWithDataFast("All");                  my ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector) = $self->RunAllStudiesWithDataFast("All");
# Line 1367  Line 1505 
1505  sub CreateSingleGenomeReactionList {  sub CreateSingleGenomeReactionList {
1506          my ($self,$RunGapFilling) = @_;          my ($self,$RunGapFilling) = @_;
1507    
1508            #Creating directory
1509            if ($self->owner() ne "master" && !-d $self->figmodel()->config("organism directory")->[0].$self->owner()."/") {
1510                    system("mkdir ".$self->figmodel()->config("organism directory")->[0].$self->owner()."/");
1511            } elsif ($self->owner() eq "master" && !-d $self->figmodel()->config("organism directory")->[0].$self->genome()."/") {
1512                    system("mkdir ".$self->figmodel()->config("organism directory")->[0].$self->genome()."/");
1513            }
1514            if ($self->owner() ne "master" && !-d $self->figmodel()->config("organism directory")->[0].$self->owner()."/".$self->genome()."/") {
1515                    system("mkdir ".$self->figmodel()->config("organism directory")->[0].$self->owner()."/".$self->genome()."/");
1516            }
1517    
1518          #Getting genome stats          #Getting genome stats
1519          my $genomestats = $self->figmodel()->get_genome_stats($self->genome());          my $genomestats = $self->figmodel()->get_genome_stats($self->genome());
1520          my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome());          my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome());
# Line 1601  Line 1749 
1749          }          }
1750    
1751          #Checking if a biomass reaction already exists          #Checking if a biomass reaction already exists
1752          my $BiomassReactionRow = $self->figmodel()->database()->get_row_by_key("BIOMASS TABLE",$self->id(),"MODELS");          my $BiomassReactionRow = $self->BuildSpecificBiomassReaction();
         if (!defined($BiomassReactionRow)) {  
                 $BiomassReactionRow = $self->BuildSpecificBiomassReaction();  
1753                  if (!defined($BiomassReactionRow)) {                  if (!defined($BiomassReactionRow)) {
1754                          $self->set_status(-2,"Preliminary reconstruction failed: could not generate biomass reaction");                          $self->set_status(-2,"Preliminary reconstruction failed: could not generate biomass reaction");
1755                          $self->figmodel()->error_message("FIGMODEL:CreateModelReactionList: Could not generate biomass function for ".$self->id()."!");                          $self->figmodel()->error_message("FIGMODEL:CreateModelReactionList: Could not generate biomass function for ".$self->id()."!");
1756                          return $self->fail();                          return $self->fail();
1757                  }                  }
         }  
1758          my $ReactionList = $BiomassReactionRow->{"ESSENTIAL REACTIONS"};          my $ReactionList = $BiomassReactionRow->{"ESSENTIAL REACTIONS"};
1759          push(@{$ReactionList},$BiomassReactionRow->{DATABASE}->[0]);          push(@{$ReactionList},$BiomassReactionRow->{DATABASE}->[0]);
1760    
# Line 1702  Line 1847 
1847          $self->figmodel()->database()->save_table($NewModelTable);          $self->figmodel()->database()->save_table($NewModelTable);
1848          $self->{_reaction_data} = $NewModelTable;          $self->{_reaction_data} = $NewModelTable;
1849          #Clearing the previous model from the cache          #Clearing the previous model from the cache
1850          $self->figmodel()->ClearDBModel($self->id(),1);          $self->figmodel()->database()->ClearDBModel($self->id(),1);
1851          #Updating the model stats table          #Updating the model stats table
1852          $self->update_stats_for_build();          $self->update_stats_for_build();
1853          $self->PrintSBMLFile();          $self->PrintSBMLFile();
# Line 1823  Line 1968 
1968          }          }
1969    
1970          #Clearing the previous model from the cache          #Clearing the previous model from the cache
1971          $self->figmodel()->ClearDBModel($self->id(),1);          $self->figmodel()->database()->ClearDBModel($self->id(),1);
1972          $ModelTable->save();          $ModelTable->save();
1973    
1974          return $self->success();          return $self->success();
# Line 1839  Line 1984 
1984  sub ArchiveModel {  sub ArchiveModel {
1985          my ($self) = @_;          my ($self) = @_;
1986    
         #Making sure the model exists  
         if (!defined($self->stats())) {  
                 $self->figmodel()->error_message("FIGMODEL:ArchiveModel: Could not find specified ".$self->id()." in database!");  
                 return $self->fail();  
         }  
   
1987          #Checking that the model file exists          #Checking that the model file exists
1988          if (!(-e $self->filename())) {          if (!(-e $self->filename())) {
1989                  $self->figmodel()->error_message("FIGMODEL:ArchiveModel: Model file ".$self->filename()." not found!");                  $self->figmodel()->error_message("FIGMODEL:ArchiveModel: Model file ".$self->filename()." not found!");
# Line 1886  Line 2025 
2025          return $self->success();          return $self->success();
2026  }  }
2027    
2028  =head2 Analysis Functions  =head2 Analysis Functions
2029    
2030  =head3 run_microarray_analysis  =head3 run_microarray_analysis
2031  Definition:  Definition:
2032          int::status = FIGMODEL->run_microarray_analysis(string::media,string::job id,string::gene calls);          int::status = FIGMODEL->run_microarray_analysis(string::media,string::job id,string::gene calls);
2033  Description:  Description:
2034          Runs microarray analysis attempting to turn off genes that are inactive in the microarray          Runs microarray analysis attempting to turn off genes that are inactive in the microarray
2035  =cut  =cut
2036  sub run_microarray_analysis {  sub run_microarray_analysis {
2037          my ($self,$media,$jobid,$index,$genecall) = @_;          my ($self,$media,$label,$index,$genecall) = @_;
2038          $genecall =~ s/_/:/g;          $genecall =~ s/_/:/g;
2039          $genecall =~ s/\//;/g;          $genecall =~ s/\//;/g;
2040          #print "\n\n".$genecall."\n\n";          my $uniqueFilename = $self->figmodel()->filename();
2041          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());          my $command = $self->figmodel()->GenerateMFAToolkitCommandLineCall($uniqueFilename,$self->id(),$media,["ProductionMFA","ShewenellaExperiment"],{"Microarray assertions" => $label.";".$index.";".$genecall,"MFASolver" => "CPLEX","Network output location" => "/scratch/"},"MicroarrayAnalysis-".$uniqueFilename.".txt",undef,$self->selected_version());
2042          #print $command."\n";          system($command);
2043          system($command);          my $filename = $self->figmodel()->config("MFAToolkit output directory")->[0].$uniqueFilename."/MicroarrayOutput-".$index.".txt";
2044          #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\"");      if (-e $filename) {
2045            my $output = $self->figmodel()->database()->load_single_column_file($filename);
2046            if (defined($output->[0])) {
2047                    my @array = split(/;/,$output->[0]);
2048                    $self->figmodel()->clearing_output($uniqueFilename,"MicroarrayAnalysis-".$uniqueFilename.".txt");
2049                    return ($array[0],$array[1],$array[8].":".$array[2],$array[9].":".$array[3],$array[10].":".$array[4],$array[11].":".$array[5],$array[12].":".$array[6],$array[13].":".$array[7]);
2050            }
2051            print STDERR $filename." is empty!";
2052        }
2053        print STDERR $filename." not found!";
2054        $self->figmodel()->clearing_output($uniqueFilename,"MicroarrayAnalysis-".$uniqueFilename.".txt");
2055    
2056            return undef;
2057    }
2058    
2059    =head3 find_minimal_pathways
2060    Definition:
2061            int::status = FIGMODEL->find_minimal_pathways(string::media,string::objective);
2062    Description:
2063            Runs microarray analysis attempting to turn off genes that are inactive in the microarray
2064    =cut
2065    sub find_minimal_pathways {
2066            my ($self,$media,$objective,$solutionnum,$AllReversible,$additionalexchange) = @_;
2067    
2068            #Setting default media
2069            if (!defined($media)) {
2070                    $media = "Complete";
2071            }
2072    
2073            #Setting default solution number
2074            if (!defined($solutionnum)) {
2075                    $solutionnum = "5";
2076            }
2077    
2078            #Setting additional exchange fluxes
2079            if (!defined($additionalexchange) || length($additionalexchange) == 0) {
2080                    if ($self->id() eq "iAF1260") {
2081                            $additionalexchange = "cpd03422[c]:-100:100;cpd01997[c]:-100:100;cpd11416[c]:-100:0;cpd15378[c]:-100:0;cpd15486[c]:-100:0";
2082                    } else {
2083                            $additionalexchange = $self->figmodel()->config("default exchange fluxes")->[0];
2084                    }
2085            }
2086    
2087            #Translating objective
2088            my $objectivestring;
2089            if ($objective eq "ALL") {
2090                    #Getting the list of universal building blocks
2091                    my $buildingblocks = $self->config("universal building blocks");
2092                    my @objectives = keys(%{$buildingblocks});
2093                    #Getting the nonuniversal building blocks
2094                    my $otherbuildingblocks = $self->config("nonuniversal building blocks");
2095                    my @array = keys(%{$otherbuildingblocks});
2096                    if (defined($self->get_biomass()) && defined($self->figmodel()->get_reaction($self->get_biomass()->{"LOAD"}->[0]))) {
2097                            my $equation = $self->figmodel()->get_reaction($self->get_biomass()->{"LOAD"}->[0])->{"EQUATION"}->[0];
2098                            if (defined($equation)) {
2099                                    for (my $i=0; $i < @array; $i++) {
2100                                            if (CORE::index($equation,$array[$i]) > 0) {
2101                                                    push(@objectives,$array[$i]);
2102                                            }
2103                                    }
2104                            }
2105                    }
2106                    for (my $i=0; $i < @objectives; $i++) {
2107                            $self->find_minimal_pathways($media,$objectives[$i]);
2108                    }
2109                    return;
2110            } elsif ($objective eq "ENERGY") {
2111                    $objectivestring = "MAX;FLUX;rxn00062;c;1";
2112            } elsif ($objective =~ m/cpd\d\d\d\d\d/) {
2113                    if ($objective =~ m/\[(\w)\]/) {
2114                            $objectivestring = "MIN;DRAIN_FLUX;".$objective.";".$1.";1";
2115                            $additionalexchange .= ";".$objective."[".$1."]:-100:0";
2116                    } else {
2117                            $objectivestring = "MIN;DRAIN_FLUX;".$objective.";c;1";
2118                            $additionalexchange .= ";".$objective."[c]:-100:0";
2119                    }
2120            } elsif ($objective =~ m/(rxn\d\d\d\d\d)/) {
2121                    my ($Reactants,$Products) = $self->figmodel()->GetReactionSubstrateData($objective);
2122                    for (my $i=0; $i < @{$Products};$i++) {
2123                            my $temp = $Products->[$i]->{"DATABASE"}->[0];
2124                            if ($additionalexchange !~ m/$temp/) {
2125                                    #$additionalexchange .= ";".$temp."[c]:-100:0";
2126                            }
2127                    }
2128                    for (my $i=0; $i < @{$Reactants};$i++) {
2129                            print $Reactants->[$i]->{"DATABASE"}->[0]." started\n";
2130                            $self->find_minimal_pathways($media,$Reactants->[$i]->{"DATABASE"}->[0],$additionalexchange);
2131                            print $Reactants->[$i]->{"DATABASE"}->[0]." done\n";
2132                    }
2133                    return;
2134            }
2135    
2136            #Adding additional drains
2137            if (($objective eq "cpd15665" || $objective eq "cpd15667" || $objective eq "cpd15668" || $objective eq "cpd15669") && $additionalexchange !~ m/cpd15666/) {
2138                    $additionalexchange .= ";cpd15666[c]:0:100";
2139            } elsif ($objective eq "cpd11493" && $additionalexchange !~ m/cpd12370/) {
2140                    $additionalexchange .= ";cpd12370[c]:0:100";
2141            } elsif ($objective eq "cpd00166" && $additionalexchange !~ m/cpd01997/) {
2142                    $additionalexchange .= ";cpd01997[c]:0:100;cpd03422[c]:0:100";
2143            }
2144    
2145            #Running MFAToolkit
2146            my $filename = $self->figmodel()->filename();
2147            my $command;
2148            if (defined($AllReversible) && $AllReversible == 1) {
2149                    $command = $self->figmodel()->GenerateMFAToolkitCommandLineCall($filename,$self->id(),$media,["ProductionMFA"],{"Make all reactions reversible in MFA"=>1, "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());
2150            } else {
2151                    $command = $self->figmodel()->GenerateMFAToolkitCommandLineCall($filename,$self->id(),$media,["ProductionMFA"],{"Make all reactions reversible in MFA"=>0, "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());
2152            }
2153            system($command);
2154    
2155            #Loading problem report
2156            my $results = $self->figmodel()->LoadProblemReport($filename);
2157            #Clearing output
2158            $self->figmodel()->clearing_output($filename,"MinimalPathways-".$media."-".$self->id()."-".$objective.".txt");
2159            if (!defined($results)) {
2160                    print STDERR $objective." pathway results not found!\n";
2161                    return;
2162            }
2163    
2164            #Parsing output
2165            my @Array;
2166            my $row = $results->get_row(1);
2167            if (defined($row->{"Notes"}->[0])) {
2168                    $_ = $row->{"Notes"}->[0];
2169                    @Array = /\d+:([^\|]+)\|/g;
2170            }
2171    
2172            #Writing output to file
2173            $self->figmodel()->database()->print_array_to_file($self->directory()."MinimalPathways-".$media."-".$objective."-".$self->id()."-".$AllReversible."-".$self->selected_version().".txt",[join("|",@Array)]);
2174  }  }
2175    
2176  =head3 find_minimal_pathways  =head3 find_minimal_pathways
# Line 1911  Line 2179 
2179  Description:  Description:
2180          Runs microarray analysis attempting to turn off genes that are inactive in the microarray          Runs microarray analysis attempting to turn off genes that are inactive in the microarray
2181  =cut  =cut
2182  sub find_minimal_pathways {  sub find_minimal_pathways_two {
2183          my ($self,$media,$objective,$solutionnum) = @_;          my ($self,$media,$objective,$solutionnum,$AllReversible,$additionalexchange) = @_;
2184    
2185          #Setting default media          #Setting default media
2186          if (!defined($media)) {          if (!defined($media)) {
# Line 1924  Line 2192 
2192                  $solutionnum = "5";                  $solutionnum = "5";
2193          }          }
2194    
2195            #Setting additional exchange fluxes
2196            if (!defined($additionalexchange) || length($additionalexchange) == 0) {
2197                    if ($self->id() eq "iAF1260") {
2198                            $additionalexchange = "cpd03422[c]:-100:100;cpd01997[c]:-100:100;cpd11416[c]:-100:0;cpd15378[c]:-100:0;cpd15486[c]:-100:0";
2199                    } else {
2200                            $additionalexchange = $self->figmodel()->config("default exchange fluxes")->[0];
2201                    }
2202            }
2203    
2204          #Translating objective          #Translating objective
2205          my $objectivestring;          my $objectivestring;
2206          if ($objective eq "ALL") {          if ($objective eq "ALL") {
# Line 1944  Line 2221 
2221                          }                          }
2222                  }                  }
2223                  for (my $i=0; $i < @objectives; $i++) {                  for (my $i=0; $i < @objectives; $i++) {
2224                          $self->find_minimal_pathways($media,$objectives[$i])                          $self->find_minimal_pathways($media,$objectives[$i]);
2225                  }                  }
2226                    return;
2227          } elsif ($objective eq "ENERGY") {          } elsif ($objective eq "ENERGY") {
2228                  $objectivestring = "MAX;FLUX;rxn00062;c;1";                  $objectivestring = "MAX;FLUX;rxn00062;c;1";
2229          } elsif ($objective =~ m/cpd\d\d\d\d\d/) {          } elsif ($objective =~ m/cpd\d\d\d\d\d/) {
2230                    if ($objective =~ m/\[(\w)\]/) {
2231                            $objectivestring = "MIN;DRAIN_FLUX;".$objective.";".$1.";1";
2232                            $additionalexchange .= ";".$objective."[".$1."]:-100:0";
2233                    } else {
2234                  $objectivestring = "MIN;DRAIN_FLUX;".$objective.";c;1";                  $objectivestring = "MIN;DRAIN_FLUX;".$objective.";c;1";
2235                            $additionalexchange .= ";".$objective."[c]:-100:0";
2236                    }
2237            } elsif ($objective =~ m/(rxn\d\d\d\d\d)/) {
2238                    my ($Reactants,$Products) = $self->figmodel()->GetReactionSubstrateData($objective);
2239                    for (my $i=0; $i < @{$Products};$i++) {
2240                            my $temp = $Products->[$i]->{"DATABASE"}->[0];
2241                            if ($additionalexchange !~ m/$temp/) {
2242                                    #$additionalexchange .= ";".$temp."[c]:-100:0";
2243                            }
2244                    }
2245                    for (my $i=0; $i < @{$Reactants};$i++) {
2246                            print $Reactants->[$i]->{"DATABASE"}->[0]." started\n";
2247                            $self->find_minimal_pathways($media,$Reactants->[$i]->{"DATABASE"}->[0],$additionalexchange);
2248                            print $Reactants->[$i]->{"DATABASE"}->[0]." done\n";
2249                    }
2250                    return;
2251          }          }
2252    
2253          #Setting additional exchange fluxes          #Adding additional drains
2254          my $additionalexchange = $self->config("default exchange fluxes")->[0];          if (($objective eq "cpd15665" || $objective eq "cpd15667" || $objective eq "cpd15668" || $objective eq "cpd15669") && $additionalexchange !~ m/cpd15666/) {
         if ($objective eq "cpd15665" || $objective eq "cpd15667" || $objective eq "cpd15668" || $objective eq "cpd15669") {  
2255                  $additionalexchange .= ";cpd15666[c]:0:100";                  $additionalexchange .= ";cpd15666[c]:0:100";
2256          } elsif ($objective eq "cpd11493") {          } elsif ($objective eq "cpd11493" && $additionalexchange !~ m/cpd12370/) {
2257                  $additionalexchange .= ";cpd12370[c]:0:100";                  $additionalexchange .= ";cpd12370[c]:0:100";
2258          } elsif ($objective eq "cpd00166") {          } elsif ($objective eq "cpd00166" && $additionalexchange !~ m/cpd01997/) {
2259                  $additionalexchange .= ";cpd01997[c]:0:100;cpd03422[c]:0:100";                  $additionalexchange .= ";cpd01997[c]:0:100;cpd03422[c]:0:100";
2260          }          }
2261    
2262          #Running MFAToolkit          #Running MFAToolkit
2263          my $filename = $self->figmodel()->filename();          my $filename = $self->figmodel()->filename();
2264          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());          my $command;
2265            if (defined($AllReversible) && $AllReversible == 1) {
2266                    $command = $self->figmodel()->GenerateMFAToolkitCommandLineCall($filename,$self->id(),$media,["ProductionMFA"],{"use simple variable and constraint names"=>1,"Make all reactions reversible in MFA"=>1, "Recursive MILP solution limit" => $solutionnum,"Generate pathways to objective" => 1,"MFASolver" => "SCIP","objective" => $objectivestring,"exchange species" => $additionalexchange},"MinimalPathways-".$media."-".$self->id().$self->selected_version().".txt",undef,$self->selected_version());
2267            } else {
2268                    $command = $self->figmodel()->GenerateMFAToolkitCommandLineCall($filename,$self->id(),$media,["ProductionMFA"],{"use simple variable and constraint names"=>1,"Make all reactions reversible in MFA"=>0, "Recursive MILP solution limit" => $solutionnum,"Generate pathways to objective" => 1,"MFASolver" => "SCIP","objective" => $objectivestring,"exchange species" => $additionalexchange},"MinimalPathways-".$media."-".$self->id().$self->selected_version().".txt",undef,$self->selected_version());
2269            }
2270            print $command."\n";
2271          system($command);          system($command);
2272    
2273          #Loading problem report          #Loading problem report
# Line 1972  Line 2275 
2275          #Clearing output          #Clearing output
2276          $self->figmodel()->clearing_output($filename,"MinimalPathways-".$media."-".$self->id()."-".$objective.".txt");          $self->figmodel()->clearing_output($filename,"MinimalPathways-".$media."-".$self->id()."-".$objective.".txt");
2277          if (!defined($results)) {          if (!defined($results)) {
2278                    print STDERR $objective." pathway results not found!\n";
2279                  return;                  return;
2280          }          }
2281    
# Line 1983  Line 2287 
2287                  @Array = /\d+:([^\|]+)\|/g;                  @Array = /\d+:([^\|]+)\|/g;
2288          }          }
2289    
2290          # Storing data in a figmodel table          #Writing output to file
2291          my $TableObject;          $self->figmodel()->database()->print_array_to_file($self->directory()."MinimalPathways-".$media."-".$objective."-".$self->id()."-".$AllReversible."-".$self->selected_version().".txt",[join("|",@Array)]);
2292          if (-e $self->directory()."MinimalPathways-".$media."-".$self->id().$self->selected_version().".txt") {  }
2293                  $TableObject->load_table($self->directory()."MinimalPathways-".$media."-".$self->id().$self->selected_version().".txt",";","|",0,["OBJECTIVE"]);  
2294          } else {  sub combine_minimal_pathways {
2295                  $TableObject = FIGMODELTable->new(["OBJECTIVE","REACTIONS"],$self->directory()."MinimalPathways-".$media."-".$self->id().$self->selected_version().".txt",["OBJECTIVE"],";","|",undef);          my ($self) = @_;
2296          }  
2297          my $tablerow = $TableObject->get_row_by_key($objective,"OBJECTIVE",1);          my $tbl;
2298          push(@{$tablerow->{"REACTIONS"}},@Array);          if (-e $self->directory()."MinimalPathwayTable-".$self->id().$self->selected_version().".tbl") {
2299          $TableObject->save();                  $tbl = FIGMODELTable::load_table($self->directory()."MinimalPathwayTable-".$self->id().$self->selected_version().".tbl",";","|",0,["Objective","Media","Reversible"]);
2300            } else {
2301                    $tbl = FIGMODELTable->new(["Objective","Media","Reactions","Reversible","Shortest path","Number of essentials","Essentials","Length"],$self->directory()."MinimalPathwayTable-".$self->id().$self->selected_version().".tbl",["Objective","Media","Reversible"],";","|");
2302            }
2303            my @files = glob($self->directory()."MinimalPathways-*");
2304            for (my $i=0; $i < @files;$i++) {
2305                    if ($files[$i] =~ m/MinimalPathways\-(\S+)\-(cpd\d\d\d\d\d)\-(\w+)\-(\d)\-/ || $files[$i] =~ m/MinimalPathways\-(\S+)\-(ENERGY)\-(\w+)\-(\d)\-/) {
2306                            my $reactions = $self->figmodel()->database()->load_single_column_file($files[$i],"");
2307                            if (defined($reactions) && @{$reactions} > 0 && length($reactions->[0]) > 0) {
2308                                    my $newrow = {"Objective"=>[$2],"Media"=>[$1],"Reversible"=>[$4]};
2309                                    my $row = $tbl->get_table_by_key($newrow->{"Objective"}->[0],"Objective")->get_table_by_key($newrow->{"Media"}->[0],"Media")->get_row_by_key($newrow->{"Reversible"}->[0],"Reversible");
2310                                    if (!defined($row)) {
2311                                            $row = $tbl->add_row($newrow);
2312                                    }
2313                                    $row->{Reactions} = $self->figmodel()->database()->load_single_column_file($files[$i],"");
2314                                    delete($row->{"Shortest path"});
2315                                    delete($row->{"Number of essentials"});
2316                                    delete($row->{"Essentials"});
2317                                    delete($row->{"Length"});
2318                                    for (my $j=0; $j < @{$row->{Reactions}}; $j++) {
2319                                            my @array = split(/,/,$row->{Reactions}->[$j]);
2320                                            $row->{"Length"}->[$j] = @array;
2321                                            if (!defined($row->{"Shortest path"}->[0]) || $row->{"Length"}->[$j] < $row->{"Shortest path"}->[0]) {
2322                                                    $row->{"Shortest path"}->[0] = $row->{"Length"}->[$j];
2323                                            }
2324                                            $row->{"Number of essentials"}->[0] = 0;
2325                                            for (my $k=0; $k < @array;$k++) {
2326                                                    if ($array[$k] =~ m/(rxn\d\d\d\d\d)/) {
2327                                                            my $class = $self->get_reaction_class($1,1);
2328                                                            my $temp = $row->{Media}->[0].":Essential";
2329                                                            if ($class =~ m/$temp/) {
2330                                                                    $row->{"Number of essentials"}->[$j]++;
2331                                                                    if (!defined($row->{"Essentials"}->[$j]) && length($row->{"Essentials"}->[$j]) > 0) {
2332                                                                            $row->{"Essentials"}->[$j] = $array[$k];
2333                                                                    } else {
2334                                                                            $row->{"Essentials"}->[$j] .= ",".$array[$k];
2335                                                                    }
2336                                                            }
2337                                                    }
2338                                            }
2339                                    }
2340                            }
2341                    }
2342            }
2343            $tbl->save();
2344  }  }
2345    
2346  =head3 calculate_growth  =head3 calculate_growth
# Line 2036  Line 2384 
2384          7.) Dead: these reactions are disconnected from the network.          7.) Dead: these reactions are disconnected from the network.
2385  =cut  =cut
2386  sub classify_model_reactions {  sub classify_model_reactions {
2387          my ($self,$Media) = @_;          my ($self,$Media,$SaveChanges) = @_;
2388    
2389          #Getting unique file for printing model output          #Getting unique file for printing model output
2390          my $UniqueFilename = $self->figmodel()->filename();          my $UniqueFilename = $self->figmodel()->filename();
# Line 2123  Line 2471 
2471                                  $CpdRow->{MEDIA}->[$index] = $Media;                                  $CpdRow->{MEDIA}->[$index] = $Media;
2472                          }                          }
2473                  }                  }
2474                    if (!defined($SaveChanges) || $SaveChanges == 1) {
2475                  $cpdclasstable->save();                  $cpdclasstable->save();
2476          }          }
2477            }
2478          if (defined($ReactionTB)) {          if (defined($ReactionTB)) {
2479                  for (my $i=0; $i < $ReactionTB->size(); $i++) {                  for (my $i=0; $i < $ReactionTB->size(); $i++) {
2480                          my $Row = $ReactionTB->get_row($i);                          my $Row = $ReactionTB->get_row($i);
# Line 2179  Line 2529 
2529                                  $RxnRow->{MEDIA}->[$index] = $Media;                                  $RxnRow->{MEDIA}->[$index] = $Media;
2530                          }                          }
2531                  }                  }
2532                    if (!defined($SaveChanges) || $SaveChanges == 1) {
2533                  $rxnclasstable->save();                  $rxnclasstable->save();
2534          }          }
2535            }
2536          return ($rxnclasstable,$cpdclasstable);          return ($rxnclasstable,$cpdclasstable);
2537  }  }
2538    
# Line 2212  Line 2564 
2564    
2565          #Running simulations          #Running simulations
2566          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");          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");
         print "simulation complete for ".$self->id().$self->selected_version()."\n";  
2567          #Parsing the results          #Parsing the results
2568          my $Results = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt",";","\\|",0,undef);          my $Results = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt",";","\\|",0,undef);
2569          if (!defined($Results)) {          if (!defined($Results)) {
# Line 2619  Line 2970 
2970                                          if (length($GapFillingRunSpecs) > 0) {                                          if (length($GapFillingRunSpecs) > 0) {
2971                                                  $GapFillingRunSpecs .= ";";                                                  $GapFillingRunSpecs .= ";";
2972                                          }                                          }
2973                                          $GapFillingRunSpecs .= $HeadingDataArray[2].":".$HeadingDataArray[1].".txt:".$HeadingDataArray[3];                                          $GapFillingRunSpecs .= $HeadingDataArray[2].":".$HeadingDataArray[1].":".$HeadingDataArray[3];
2974                                  } else {                                  } else {
2975                                          $SolutionExistedCount++;                                          $SolutionExistedCount++;
2976                                  }                                  }
# Line 2644  Line 2995 
2995          my $SolutionsFound = 0;          my $SolutionsFound = 0;
2996          my $GapFillingArray;          my $GapFillingArray;
2997          push(@{$GapFillingArray},split(/;/,$GapFillingRunSpecs));          push(@{$GapFillingArray},split(/;/,$GapFillingRunSpecs));
2998          $self->datagapfill($GapFillingArray);          my $GapFillingResults = $self->datagapfill($GapFillingArray,"GFS");
2999            if (defined($GapFillingResults)) {
3000                    $SolutionsFound = 1;
3001            }
3002    
3003          if (defined($RescuedPreviousResults) && @{$RescuedPreviousResults} > 0) {          if (defined($RescuedPreviousResults) && @{$RescuedPreviousResults} > 0) {
3004                  #Printing previous solutions to GFS file                  #Printing previous solutions to GFS file
# Line 2667  Line 3021 
3021          return $self->success();          return $self->success();
3022  }  }
3023    
3024    =head3 SolutionReconciliation
3025    Definition:
3026            FIGMODELmodel->SolutionReconciliation();
3027    Description:
3028            This is a wrapper for running the solution reconciliation algorithm on any model in the database.
3029            The algorithm performs a reconciliation of any gap filling solutions to identify the combination of solutions that results in the optimal model.
3030            This function prints out one output file in the Model directory: ReconciliationOutput.txt: this is a summary of the results of the reconciliation analysis
3031    =cut
3032    
3033    sub SolutionReconciliation {
3034            my ($self,$GapFill,$Stage) = @_;
3035    
3036            #Setting the output filenames
3037            my $OutputFilename;
3038            my $OutputFilenameTwo;
3039            if ($GapFill == 1) {
3040                    $OutputFilename = $self->directory().$self->id().$self->selected_version()."-GFReconciliation.txt";
3041                    $OutputFilenameTwo = $self->directory().$self->id().$self->selected_version()."-GFSRS.txt";
3042            } else {
3043                    $OutputFilename = $self->directory().$self->id().$self->selected_version()."-GGReconciliation.txt";
3044                    $OutputFilenameTwo = $self->directory().$self->id().$self->selected_version()."-GGSRS.txt";
3045            }
3046    
3047            #In stage one, we run the reconciliation and create a test file to check combined solution performance
3048            if (!defined($Stage) || $Stage == 1) {
3049                    my $GrowMatchTable = $self->figmodel()->database()->LockDBTable("GROWMATCH TABLE");
3050                    my $Row = $GrowMatchTable->get_row_by_key($self->genome(),"ORGANISM",1);
3051                    $Row->{"GF RECONCILATION TIMING"}->[0] = time()."-";
3052                    $GrowMatchTable->save();
3053                    $self->figmodel()->database()->UnlockDBTable("GROWMATCH TABLE");
3054    
3055                    #Getting a unique filename
3056                    my $UniqueFilename = $self->figmodel()->filename();
3057    
3058                    #Copying over the necessary files
3059                    if ($GapFill == 1) {
3060                            if (!-e $self->directory().$self->id().$self->selected_version()."-GFEM.txt") {
3061                                    print STDERR "FIGMODEL:SolutionReconciliation:".$self->directory().$self->id().$self->selected_version()."-GFEM.txt file not found. Could not reconcile!";
3062                                    return 0;
3063                            }
3064                            if (!-e $self->directory().$self->id().$self->selected_version()."-OPEM.txt") {
3065                                    print STDERR "FIGMODEL:SolutionReconciliation:".$self->directory().$self->id().$self->selected_version()."-OPEM.txt file not found. Could not reconcile!";
3066                                    return 0;
3067                            }
3068                            system("cp ".$self->directory().$self->id().$self->selected_version()."-GFEM.txt ".$self->figmodel()->config("MFAToolkit input files")->[0].$UniqueFilename."-GFEM.txt");
3069                            system("cp ".$self->directory().$self->id().$self->selected_version()."-OPEM.txt ".$self->figmodel()->config("MFAToolkit input files")->[0].$UniqueFilename."-OPEM.txt");
3070                            #Backing up and deleting the existing reconciliation file
3071                            if (-e $OutputFilename) {
3072                                    system("cp ".$OutputFilename." ".$self->directory().$self->id().$self->selected_version()."-OldGFReconciliation.txt");
3073                                    unlink($OutputFilename);
3074                            }
3075                    } else {
3076                            if (!-e $self->directory().$self->id().$self->selected_version()."-GGEM.txt") {
3077                                    print STDERR "FIGMODEL:SolutionReconciliation:".$self->directory().$self->id().$self->selected_version()."-GGEM.txt file not found. Could not reconcile!";
3078                                    return 0;
3079                            }
3080                            if (!-e $self->directory().$self->id().$self->selected_version()."-GGOPEM.txt") {
3081                                    print STDERR "FIGMODEL:SolutionReconciliation:".$self->directory().$self->id().$self->selected_version()."-GGOPEM.txt file not found. Could not reconcile!";
3082                                    return 0;
3083                            }
3084                            system("cp ".$self->directory().$self->id().$self->selected_version()."-GGEM.txt ".$self->figmodel()->config("MFAToolkit input files")->[0].$UniqueFilename."-GGEM.txt");
3085                            system("cp ".$self->directory().$self->id().$self->selected_version()."-GGOPEM.txt ".$self->figmodel()->config("MFAToolkit input files")->[0].$UniqueFilename."-OPEM.txt");
3086                            #Backing up and deleting the existing reconciliation file
3087                            if (-e $OutputFilename) {
3088                                    system("cp ".$OutputFilename." ".$self->directory().$self->id().$self->selected_version()."-OldGGReconciliation.txt");
3089                                    unlink($OutputFilename);
3090                            }
3091                    }
3092    
3093                    #Running the reconciliation
3094                    system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),"NONE",["SolutionReconciliation"],{"Solution data for model optimization" => $UniqueFilename},"Reconciliation".$UniqueFilename.".log",undef,$self->selected_version()));
3095                    $GrowMatchTable = $self->figmodel()->database()->LockDBTable("GROWMATCH TABLE");
3096                    $Row = $GrowMatchTable->get_row_by_key($self->genome(),"ORGANISM",1);
3097                    $Row->{"GF RECONCILATION TIMING"}->[0] .= time();
3098                    $GrowMatchTable->save();
3099                    $self->figmodel()->database()->UnlockDBTable("GROWMATCH TABLE");
3100    
3101                    #Loading the problem report from the reconciliation run
3102                    my $ReconciliatonOutput = $self->figmodel()->LoadProblemReport($UniqueFilename);
3103                    print $UniqueFilename."\n";
3104                    #Clearing output files
3105                    $self->figmodel()->clearing_output($UniqueFilename,"Reconciliation".$UniqueFilename.".log");
3106                    $ReconciliatonOutput->save("/home/chenry/Test.txt");
3107    
3108                    #Checking the a problem report was found and was loaded
3109                    if (!defined($ReconciliatonOutput) || $ReconciliatonOutput->size() < 1 || !defined($ReconciliatonOutput->get_row(0)->{"Notes"}->[0])) {
3110                            print STDERR "FIGMODEL:SolutionReconciliation: MFAToolkit output from SolutionReconciliation of ".$self->id()." not found!\n\n";
3111                            return 0;
3112                    }
3113    
3114                    #Processing the solutions
3115                    my $SolutionCount = 0;
3116                    my $ReactionSetHash;
3117                    my $SingleReactionHash;
3118                    my $ReactionDataHash;
3119                    for (my $n=0; $n < $ReconciliatonOutput->size(); $n++) {
3120                            if (defined($ReconciliatonOutput->get_row($n)->{"Notes"}->[0]) && $ReconciliatonOutput->get_row($n)->{"Notes"}->[0] =~ m/^Recursive\sMILP\s([^;]+)/) {
3121                                    #Breaking up the solution into reaction sets
3122                                    my @ReactionSets = split(/\|/,$1);
3123                                    #Creating reaction lists for each set
3124                                    my $SolutionHash;
3125                                    for (my $i=0; $i < @ReactionSets; $i++) {
3126                                            if (length($ReactionSets[$i]) > 0) {
3127                                                    my @Alternatives = split(/:/,$ReactionSets[$i]);
3128                                                    for (my $j=1; $j < @Alternatives; $j++) {
3129                                                            if (length($Alternatives[$j]) > 0) {
3130                                                                    push(@{$SolutionHash->{$Alternatives[$j]}},$Alternatives[0]);
3131                                                            }
3132                                                    }
3133                                                    if (@Alternatives == 1) {
3134                                                            $SingleReactionHash->{$Alternatives[0]}->{$SolutionCount} = 1;
3135                                                            if (!defined($SingleReactionHash->{$Alternatives[0]}->{"COUNT"})) {
3136                                                                    $SingleReactionHash->{$Alternatives[0]}->{"COUNT"} = 0;
3137                                                            }
3138                                                            $SingleReactionHash->{$Alternatives[0]}->{"COUNT"}++;
3139                                                    }
3140                                            }
3141                                    }
3142                                    #Identifying reactions sets and storing the sets in the reactions set hash
3143                                    foreach my $Solution (keys(%{$SolutionHash})) {
3144                                            my $SetKey = join(",",sort(@{$SolutionHash->{$Solution}}));
3145                                            if (!defined($ReactionSetHash->{$SetKey}->{$SetKey}->{$SolutionCount})) {
3146                                                    $ReactionSetHash->{$SetKey}->{$SetKey}->{$SolutionCount} = 1;
3147                                                    if (!defined($ReactionSetHash->{$SetKey}->{$SetKey}->{"COUNT"})) {
3148                                                            $ReactionSetHash->{$SetKey}->{$SetKey}->{"COUNT"} = 0;
3149                                                    }
3150                                                    $ReactionSetHash->{$SetKey}->{$SetKey}->{"COUNT"}++;
3151                                            }
3152                                            $ReactionSetHash->{$SetKey}->{$Solution}->{$SolutionCount} = 1;
3153                                            if (!defined($ReactionSetHash->{$SetKey}->{$Solution}->{"COUNT"})) {
3154                                                    $ReactionSetHash->{$SetKey}->{$Solution}->{"COUNT"} = 0;
3155                                            }
3156                                            $ReactionSetHash->{$SetKey}->{$Solution}->{"COUNT"}++;
3157                                    }
3158                                    $SolutionCount++;
3159                            }
3160                    }
3161    
3162                    #Handling the scenario where no solutions were found
3163                    if ($SolutionCount == 0) {
3164                            print STDERR "FIGMODEL:SolutionReconciliation: Reconciliation unsuccessful. No solution found.\n\n";
3165                            return 0;
3166                    }
3167    
3168                    #Printing results without solution performance figures. Also printing solution test file
3169                    open (RECONCILIATION, ">$OutputFilename");
3170                    #Printing the file heading
3171                    print RECONCILIATION "DATABASE;DEFINITION;REVERSIBLITY;DELTAG;DIRECTION;NUMBER OF SOLUTIONS";
3172                    for (my $i=0; $i < $SolutionCount; $i++) {
3173                            print RECONCILIATION ";Solution ".$i;
3174                    }
3175                    print RECONCILIATION "\n";
3176                    #Printing the singlet reactions first
3177                    my $Solutions;
3178                    print RECONCILIATION "SINGLET REACTIONS\n";
3179                    my @SingletReactions = keys(%{$SingleReactionHash});
3180                    for (my $j=0; $j < $SolutionCount; $j++) {
3181                            $Solutions->[$j]->{"BASE"} = $j;
3182                    }
3183                    for (my $i=0; $i < @SingletReactions; $i++) {
3184                            my $ReactionData;
3185                            if (defined($ReactionDataHash->{$SingletReactions[$i]})) {
3186                                    $ReactionData = $ReactionDataHash->{$SingletReactions[$i]};
3187                            } else {
3188                                    my $Direction = substr($SingletReactions[$i],0,1);
3189                                    if ($Direction eq "+") {
3190                                            $Direction = "=>";
3191                                    } else {
3192                                            $Direction = "<=";
3193                                    }
3194                                    my $Reaction = substr($SingletReactions[$i],1);
3195                                    $ReactionData = FIGMODELObject->load($self->figmodel()->config("reaction directory")->[0].$Reaction,"\t");
3196                                    $ReactionData->{"DIRECTIONS"}->[0] = $Direction;
3197                                    $ReactionData->{"REACTIONS"}->[0] = $Reaction;
3198                                    if (!defined($ReactionData->{"DEFINITION"}->[0])) {
3199                                            $ReactionData->{"DEFINITION"}->[0] = "UNKNOWN";
3200                                    }
3201                                    if (!defined($ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0])) {
3202                                            $ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0] = "UNKNOWN";
3203                                    }
3204                                    if (!defined($ReactionData->{"DELTAG"}->[0])) {
3205                                            $ReactionData->{"DELTAG"}->[0] = "UNKNOWN";
3206                                    }
3207                                    $ReactionDataHash->{$SingletReactions[$i]} = $ReactionData;
3208                            }
3209                            print RECONCILIATION $ReactionData->{"REACTIONS"}->[0].";".$ReactionData->{"DEFINITION"}->[0].";".$ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0].";".$ReactionData->{"DELTAG"}->[0].";".$ReactionData->{"DIRECTIONS"}->[0].";".$SingleReactionHash->{$SingletReactions[$i]}->{"COUNT"};
3210                            for (my $j=0; $j < $SolutionCount; $j++) {
3211                                    print RECONCILIATION ";";
3212                                    if (defined($SingleReactionHash->{$SingletReactions[$i]}->{$j})) {
3213                                            $Solutions->[$j]->{$SingletReactions[$i]} = 1;
3214                                            $Solutions->[$j]->{"BASE"} = $j;
3215                                            print RECONCILIATION "|".$j."|";
3216                                    }
3217                            }
3218                            print RECONCILIATION "\n";
3219                    }
3220                    #Printing the reaction sets with alternatives
3221                    print RECONCILIATION "Reaction sets with alternatives\n";
3222                    my @ReactionSets = keys(%{$ReactionSetHash});
3223                    foreach my $ReactionSet (@ReactionSets) {
3224                            my $NewSolutions;
3225                            my $BaseReactions;
3226                            my $AltList = [$ReactionSet];
3227                            push(@{$AltList},keys(%{$ReactionSetHash->{$ReactionSet}}));
3228                            for (my $j=0; $j < @{$AltList}; $j++) {
3229                                    my $CurrentNewSolutions;
3230                                    my $Index;
3231                                    if ($j == 0) {
3232                                            print RECONCILIATION "NEW SET\n";
3233                                    } elsif ($AltList->[$j] ne $ReactionSet) {
3234                                            print RECONCILIATION "ALTERNATIVE SET\n";
3235                                            #For each base solution in which this set is represented, we copy the base solution to the new solution
3236                                            my $NewSolutionCount = 0;
3237                                            for (my $k=0; $k < $SolutionCount; $k++) {
3238                                                    if (defined($ReactionSetHash->{$ReactionSet}->{$AltList->[$j]}->{$k})) {
3239                                                            if (defined($Solutions)) {
3240                                                                    $Index->{$k} = @{$Solutions} + $NewSolutionCount;
3241                                                            } else {
3242                                                                    $Index->{$k} = $NewSolutionCount;
3243                                                            }
3244                                                            if (defined($NewSolutions) && @{$NewSolutions} > 0) {
3245                                                                    $Index->{$k} += @{$NewSolutions};
3246                                                            }
3247                                                            $CurrentNewSolutions->[$NewSolutionCount] = {};
3248                                                            foreach my $Reaction (keys(%{$Solutions->[$k]})) {
3249                                                                    $CurrentNewSolutions->[$NewSolutionCount]->{$Reaction} = $Solutions->[$k]->{$Reaction};
3250                                                            }
3251                                                            $NewSolutionCount++;
3252                                                    }
3253                                            }
3254                                    }
3255                                    if ($j == 0 || $AltList->[$j] ne $ReactionSet) {
3256                                            my @SingletReactions = split(/,/,$AltList->[$j]);
3257                                            for (my $i=0; $i < @SingletReactions; $i++) {
3258                                                    #Adding base reactions to base solutions and set reactions the new solutions
3259                                                    if ($j == 0) {
3260                                                            push(@{$BaseReactions},$SingletReactions[$i]);
3261                                                    } else {
3262                                                            for (my $k=0; $k < @{$CurrentNewSolutions}; $k++) {
3263                                                                    $CurrentNewSolutions->[$k]->{$SingletReactions[$i]} = 1;
3264                                                            }
3265                                                    }
3266                                                    #Getting reaction data and printing reaction in output file
3267                                                    my $ReactionData;
3268                                                    if (defined($ReactionDataHash->{$SingletReactions[$i]})) {
3269                                                            $ReactionData = $ReactionDataHash->{$SingletReactions[$i]};
3270                                                    } else {
3271                                                            my $Direction = substr($SingletReactions[$i],0,1);
3272                                                            if ($Direction eq "+") {
3273                                                                    $Direction = "=>";
3274                                                            } else {
3275                                                                    $Direction = "<=";
3276                                                            }
3277                                                            my $Reaction = substr($SingletReactions[$i],1);
3278                                                            $ReactionData = FIGMODELObject->load($self->figmodel()->config("reaction directory")->[0].$Reaction,"\t");
3279                                                            $ReactionData->{"DIRECTIONS"}->[0] = $Direction;
3280                                                            $ReactionData->{"REACTIONS"}->[0] = $Reaction;
3281                                                            if (!defined($ReactionData->{"DEFINITION"}->[0])) {
3282                                                                    $ReactionData->{"DEFINITION"}->[0] = "UNKNOWN";
3283                                                            }
3284                                                            if (!defined($ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0])) {
3285                                                                    $ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0] = "UNKNOWN";
3286                                                            }
3287                                                            if (!defined($ReactionData->{"DELTAG"}->[0])) {
3288                                                                    $ReactionData->{"DELTAG"}->[0] = "UNKNOWN";
3289                                                            }
3290                                                            $ReactionDataHash->{$SingletReactions[$i]} = $ReactionData;
3291                                                    }
3292                                                    print RECONCILIATION $ReactionData->{"REACTIONS"}->[0].";".$ReactionData->{"DEFINITION"}->[0].";".$ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0].";".$ReactionData->{"DELTAG"}->[0].";".$ReactionData->{"DIRECTIONS"}->[0].";".$ReactionSetHash->{$ReactionSet}->{$AltList->[$j]}->{"COUNT"};
3293                                                    for (my $k=0; $k < $SolutionCount; $k++) {
3294                                                            print RECONCILIATION ";";
3295                                                            if (defined($ReactionSetHash->{$ReactionSet}->{$AltList->[$j]}->{$k})) {
3296                                                                    if ($j == 0) {
3297                                                                            print RECONCILIATION "|".$k."|";
3298                                                                    } else {
3299                                                                            print RECONCILIATION "|".$Index->{$k}."|";
3300                                                                    }
3301                                                            }
3302                                                    }
3303                                                    print RECONCILIATION "\n";
3304                                            }
3305                                            #Adding the current new solutions to the new solutions array
3306                                            if (defined($CurrentNewSolutions) && @{$CurrentNewSolutions} > 0) {
3307                                                    push(@{$NewSolutions},@{$CurrentNewSolutions});
3308                                            }
3309                                    }
3310                            }
3311                            #Adding the base reactions to all existing solutions
3312                            for (my $j=0; $j < @{$Solutions}; $j++) {
3313                                    if (defined($ReactionSetHash->{$ReactionSet}->{$ReactionSet}->{$Solutions->[$j]->{"BASE"}})) {
3314                                            foreach my $SingleReaction (@{$BaseReactions}) {
3315                                                    $Solutions->[$j]->{$SingleReaction} = 1;
3316                                            }
3317                                    }
3318                            }
3319                            #Adding the new solutions to the set of existing solutions
3320                            push(@{$Solutions},@{$NewSolutions});
3321                    }
3322                    close(RECONCILIATION);
3323                    #Now printing a file that defines all of the solutions in a format the testsolutions function understands
3324                    open (RECONCILIATION, ">$OutputFilenameTwo");
3325                    print RECONCILIATION "Experiment;Solution index;Solution cost;Solution reactions\n";
3326                    for (my $i=0; $i < @{$Solutions}; $i++) {
3327                            delete($Solutions->[$i]->{"BASE"});
3328                            print RECONCILIATION "SR".$i.";".$i.";10;".join(",",keys(%{$Solutions->[$i]}))."\n";
3329                    }
3330                    close(RECONCILIATION);
3331    
3332                    $GrowMatchTable = $self->figmodel()->database()->LockDBTable("GROWMATCH TABLE");
3333                    $Row = $GrowMatchTable->get_row_by_key($self->genome(),"ORGANISM",1);
3334                    $Row->{"GF RECON TESTING TIMING"}->[0] = time()."-";
3335                    $Row->{"GF RECON SOLUTIONS"}->[0] = @{$Solutions};
3336                    $GrowMatchTable->save();
3337                    $self->figmodel()->database()->UnlockDBTable("GROWMATCH TABLE");
3338    
3339                    #Scheduling the solution testing
3340                    if ($GapFill == 1) {
3341                            system($self->figmodel()->config("scheduler executable")->[0]." \"add:testsolutions?".$self->id().$self->selected_version()."?-1?GFSR:BACK:fast:QSUB\"");
3342                    } else {
3343                            system($self->figmodel()->config("scheduler executable")->[0]." \"add:testsolutions?".$self->id().$self->selected_version()."?-1?GGSR:BACK:fast:QSUB\"");
3344                    }
3345            } else {
3346                    #Reading in the solution testing results
3347                    my $Data;
3348                    if ($GapFill == 1) {
3349                            $Data = $self->figmodel()->database()->load_single_column_file($self->directory().$self->id().$self->selected_version()."-GFSREM.txt","");
3350                    } else {
3351                            $Data = $self->figmodel()->database()->load_single_column_file($self->directory().$self->id().$self->selected_version()."-GGSREM.txt","");
3352                    }
3353    
3354                    #Reading in the preliminate reconciliation report
3355                    my $OutputData = $self->figmodel()->database()->load_single_column_file($OutputFilename,"");
3356                    #Replacing the file tags with actual performance data
3357                    my $Count = 0;
3358                    for (my $i=0; $i < @{$Data}; $i++) {
3359                            if ($Data->[$i] =~ m/^SR(\d+);.+;(\d+\/\d+);/) {
3360                                    my $Index = $1;
3361                                    my $Performance = $Index."/".$2;
3362                                    for (my $j=0; $j < @{$OutputData}; $j++) {
3363                                            $OutputData->[$j] =~ s/\|$Index\|/$Performance/g;
3364                                    }
3365                            }
3366                    }
3367                    $self->figmodel()->database()->print_array_to_file($OutputFilename,$OutputData);
3368    
3369                    my $GrowMatchTable = $self->figmodel()->database()->LockDBTable("GROWMATCH TABLE");
3370                    my $Row = $GrowMatchTable->get_row_by_key($self->genome(),"ORGANISM",1);
3371                    $Row->{"GF RECON TESTING TIMING"}->[0] .= time();
3372                    $GrowMatchTable->save();
3373                    $self->figmodel()->database()->UnlockDBTable("GROWMATCH TABLE");
3374            }
3375    
3376            return 1;
3377    }
3378    
3379  =head3 BuildSpecificBiomassReaction  =head3 BuildSpecificBiomassReaction
3380  Definition:  Definition:
3381          FIGMODELmodel->BuildSpecificBiomassReaction();          FIGMODELmodel->BuildSpecificBiomassReaction();
# Line 2679  Line 3388 
3388          my $OrganismID = $self->genome();          my $OrganismID = $self->genome();
3389          #Checking for a biomass override          #Checking for a biomass override
3390          if (defined($self->config("biomass reaction override")->{$OrganismID})) {          if (defined($self->config("biomass reaction override")->{$OrganismID})) {
3391                  $biomassrxn = $self->config("biomass reaction override")->{$OrganismID};                  my $biomassID = $self->config("biomass reaction override")->{$OrganismID};
3392                  print "Overriding biomass template and selecting ".$biomassrxn." for ".$OrganismID.".\n";                  my $tbl = $self->database()->get_table("BIOMASS",1);
3393                    $biomassrxn = $tbl->get_row_by_key($biomassID,"DATABASE");
3394                    print "Overriding biomass template and selecting ".$biomassID." for ".$OrganismID.".\n";
3395          } else {#Creating biomass reaction from the template          } else {#Creating biomass reaction from the template
3396                  #Getting the genome stats                  #Getting the genome stats
3397                  my $genomestats = $self->figmodel()->get_genome_stats($self->genome());                  my $genomestats = $self->figmodel()->get_genome_stats($self->genome());
# Line 2971  Line 3682 
3682    
3683                  #Adding the biomass equation to the biomass table                  #Adding the biomass equation to the biomass table
3684                  my $NewRow = $self->figmodel()->add_biomass_reaction($Equation,$self->id(),"Template:".$self->genome());                  my $NewRow = $self->figmodel()->add_biomass_reaction($Equation,$self->id(),"Template:".$self->genome());
3685                  $biomassrxn = $NewRow->{DATABASE}->[0];                  $biomassrxn = $NewRow;
                 print $biomassrxn."\n";  
3686          }          }
3687          print $biomassrxn."\n";          return $biomassrxn;
         my $BiomassRow = $self->figmodel()->add_model_to_biomass_reaction($biomassrxn,$self->id());  
         return $BiomassRow;  
3688  }  }
3689    
3690  =head3 PrintSBMLFile  =head3 PrintSBMLFile
# Line 3287  Line 3995 
3995          $self->figmodel()->clearing_output($UniqueFilename,"FBA-".$self->id().$self->selected_version().".lp");          $self->figmodel()->clearing_output($UniqueFilename,"FBA-".$self->id().$self->selected_version().".lp");
3996  }  }
3997    
3998    =head3 patch_model
3999    Definition:
4000            FIGMODELTable:patch results = FIGMODELmodel->patch_model(FIGMODELTable:patch table);
4001    Description:
4002    =cut
4003    sub patch_model {
4004            my ($self,$tbl) = @_;
4005    
4006            #Instantiating table
4007            my $results = FIGMODELTable->new(["Reactions","New genes","Old genes","Genes","Roles","Status"],$self->directory()."PatchResults-".$self->id().$self->selected_version().".tbl",["Reaction"],"\t",";",undef);
4008            #Getting genome annotations
4009            my $features = $self->figmodel()->database()->get_genome_feature_table($self->genome());
4010            #Gettubg reaction table
4011            my $reactions = $self->reaction_table();
4012            #Checking for patched roles
4013            for (my $i=0; $i < $tbl->size(); $i++) {
4014                    my $row = $tbl->get_row($i);
4015                    my @genes = $features->get_rows_by_key($row->{ROLE}->[0],"ROLES");
4016                    if (@genes > 0) {
4017                            for (my $j=0; $j < @{$row->{REACTIONS}};$j++) {
4018                                    my $resultrxn = $results->get_row_by_key($row->{REACTIONS}->[$j],"Reactions");
4019                                    if (!defined($resultrxn)) {
4020                                            $resultrxn = $results->add_row({"Reactions"=>[$row->{REACTIONS}->[$j]],"Roles"=>[$row->{ROLE}->[0]]});
4021                                    }
4022                                    my $rxnrow = $reactions->get_row_by_key($row->{REACTIONS}->[$j],"LOAD");
4023                                    if (defined($rxnrow) && !defined($resultrxn->{"Old genes"})) {
4024                                            $resultrxn->{"Old genes"} = $rxnrow->{"ASSOCIATED PEG"};
4025                                            if ($resultrxn->{"Old genes"}->[0] !~ m/GAP|BOF|UNIVERSAL|SPONTANEOUS/) {
4026                                                    push(@{$resultrxn->{"Genes"}},@{$resultrxn->{"Old genes"}});
4027                                            }
4028                                    }
4029                                    delete $resultrxn->{"Current gene set"};
4030                                    if (defined($resultrxn->{"Genes"})) {
4031                                            push(@{$resultrxn->{"Current gene set"}},@{$resultrxn->{"Genes"}});
4032                                    }
4033                                    for (my $k=0; $k < @genes; $k++) {
4034                                            if ($genes[$k]->{ID}->[0] =~ m/(peg\.\d+)/) {
4035                                                    my $gene = $1;
4036                                                    my $addgene = 1;
4037                                                    if (defined($resultrxn->{"Old genes"})) {
4038                                                            for (my $m=0; $m < @{$resultrxn->{"Old genes"}}; $m++) {
4039                                                                    if ($resultrxn->{"Old genes"}->[$m] =~ m/$gene/) {
4040                                                                            $addgene = 0;
4041                                                                    }
4042                                                            }
4043                                                    }
4044                                                    if ($addgene == 1) {
4045                                                            push(@{$resultrxn->{"New genes"}},$gene);
4046                                                            if ($row->{COMPLEX}->[0] ne "0" && defined($resultrxn->{"Current gene set"})) {
4047                                                                    my $added = 0;
4048                                                                    for (my $m=0; $m < @{$resultrxn->{"Current gene set"}}; $m++) {
4049                                                                            if ($row->{COMPLEX}->[0] eq "1") {
4050                                                                                    $resultrxn->{"Current gene set"}->[$m] = $resultrxn->{"Current gene set"}->[$m]."+".$gene;
4051                                                                                    $added = 1;
4052                                                                            } else {
4053                                                                                    my @geneset = split(/\+/,$resultrxn->{"Current gene set"}->[$m]);
4054                                                                                    for (my $n=0; $n < @geneset;$n++) {
4055                                                                                            if ($self->figmodel()->colocalized_genes($geneset[$n],$gene,$self->genome()) == 1) {
4056                                                                                                    $resultrxn->{"Current gene set"}->[$m] = $resultrxn->{"Current gene set"}->[$m]."+".$gene;
4057                                                                                                    $added = 1;
4058                                                                                                    last;
4059                                                                                            }
4060                                                                                    }
4061                                                                            }
4062                                                                    }
4063                                                                    if ($added == 0) {
4064                                                                            push(@{$resultrxn->{"Current gene set"}},$gene);
4065                                                                    }
4066                                                            } else {
4067                                                                    push(@{$resultrxn->{"Current gene set"}},$gene);
4068                                                            }
4069                                                    }
4070                                            }
4071                                    }
4072                                    delete $resultrxn->{"Genes"};
4073                                    push(@{$resultrxn->{"Genes"}},@{$resultrxn->{"Current gene set"}});
4074                            }
4075                    }
4076            }
4077    
4078            #Ensuring that the old model is preserved
4079            $self->ArchiveModel();
4080            #Modifing the reaction list
4081            for (my $i=0; $i < $results->size();$i++) {
4082                    my $row = $results->get_row($i);
4083                    my $rxnrow = $reactions->get_row_by_key($row->{"Reactions"}->[0],"LOAD");
4084                    if (defined($rxnrow)) {
4085                            $rxnrow->{"ASSOCIATED PEG"} = $row->{"Genes"};
4086                    } else {
4087                            $reactions->add_row({LOAD=>[$row->{"Reactions"}->[0]],DIRECTIONALITY=>[$self->figmodel()->reversibility_of_reaction($row->{"Reactions"}->[0])],COMPARTMENT=>["c"],"ASSOCIATED PEG"=>$row->{"Genes"},SUBSYSTEM=>["NONE"],CONFIDENCE=>[2],REFERENCE=>["NONE"],NOTES=>["PATCH"]});
4088                    }
4089            }
4090            $reactions->save();
4091            $results->save();
4092            $self->update_model_stats();
4093            $self->PrintModelLPFile();
4094            $self->run_default_model_predictions();
4095            #Returning results
4096            return $results;
4097    }
4098    
4099    =head3 translate_genes
4100    Definition:
4101            FIGMODELmodel->translate_genes();
4102    Description:
4103    =cut
4104    sub translate_genes {
4105            my ($self) = @_;
4106    
4107            #Loading gene translations
4108            if (!defined($self->{_gene_aliases})) {
4109                    #Loading gene aliases from feature table
4110                    my $tbl = $self->figmodel()->GetGenomeFeatureTable($self->genome());
4111                    if (defined($tbl)) {
4112                            for (my $i=0; $i < $tbl->size(); $i++) {
4113                                    my $row = $tbl->get_row($i);
4114                                    if ($row->{ID}->[0] =~ m/(peg\.\d+)/) {
4115                                            my $geneID = $1;
4116                                            for (my $j=0; $j < @{$row->{ALIASES}}; $j++) {
4117                                                    $self->{_gene_aliases}->{$row->{ALIASES}->[$j]} = $geneID;
4118                                            }
4119                                    }
4120                            }
4121                    }
4122                    #Loading additional gene aliases from the database
4123                    if (-e $self->figmodel()->config("Translation directory")->[0]."AdditionalAliases/".$self->genome().".txt") {
4124                            my $AdditionalAliases = $self->figmodel()->database()->load_multiple_column_file($self->figmodel()->config("Translation directory")->[0]."AdditionalAliases/".$self->genome().".txt","\t");
4125                            for (my $i=0; $i < @{$AdditionalAliases}; $i++) {
4126                                    $self->{_gene_aliases}->{$AdditionalAliases->[$i]->[1]} = $AdditionalAliases->[$i]->[0];
4127                            }
4128                    }
4129            }
4130    
4131            #Cycling through reactions and translating genes
4132            for (my $i=0; $i < $self->reaction_table()->size(); $i++) {
4133                    my $row = $self->reaction_table()->get_row($i);
4134                    if (defined($row->{"ASSOCIATED PEG"})) {
4135                            for (my $j=0; $j < @{$row->{"ASSOCIATED PEG"}}; $j++) {
4136                                    my $Original = $row->{"ASSOCIATED PEG"}->[$j];
4137                                    $Original =~ s/\sand\s/:/g;
4138                                    $Original =~ s/\sor\s/;/g;
4139                                    my @GeneNames = split(/[,\+\s\(\):;]/,$Original);
4140                                    foreach my $Gene (@GeneNames) {
4141                                            if (length($Gene) > 0 && defined($self->{_gene_aliases}->{$Gene})) {
4142                                                    my $Replace = $self->{_gene_aliases}->{$Gene};
4143                                                    $Original =~ s/([^\w])$Gene([^\w])/$1$Replace$2/g;
4144                                                    $Original =~ s/^$Gene([^\w])/$Replace$1/g;
4145                                                    $Original =~ s/([^\w])$Gene$/$1$Replace/g;
4146                                                    $Original =~ s/^$Gene$/$Replace/g;
4147                                            }
4148                                    }
4149                                    $Original =~ s/:/ and /g;
4150                                    $Original =~ s/;/ or /g;
4151                                    $row->{"ASSOCIATED PEG"}->[$j] = $Original;
4152                            }
4153                    }
4154            }
4155    
4156            #Archiving model and saving reaction table
4157            $self->ArchiveModel();
4158            $self->reaction_table()->save();
4159    }
4160    
4161    =head3 feature_web_data
4162    Definition:
4163            string:web output for feature/model connection = FIGMODELmodel->feature_web_data(FIGMODELfeature:feature);
4164    Description:
4165    =cut
4166    sub feature_web_data {
4167            my ($self,$feature) = @_;
4168    
4169            #First checking if the feature is in the model
4170            my $data = $self->get_feature_data($feature->{ID}->[0]);
4171            if (!defined($data)) {
4172                    return "Not in model";
4173            }
4174    
4175            my $output;
4176            #Printing predictions
4177            if (defined($data->{$self->id()."PREDICTIONS"})) {
4178                    #Parsing essentiality data
4179                    my $essentialityData;
4180                    if (defined($feature->{ESSENTIALITY})) {
4181                            for (my $i=0; $i < @{$feature->{ESSENTIALITY}};$i++) {
4182                                    my @temp = split(/:/,$feature->{ESSENTIALITY}->[$i]);
4183                                    $essentialityData->{$temp[0]} = $temp[1];
4184                            }
4185                    }
4186                    my $predictionHash;
4187                    for (my $i=0; $i < @{$data->{$self->id()."PREDICTIONS"}};$i++) {
4188                            my @temp = split(/:/,$data->{$self->id()."PREDICTIONS"}->[$i]);
4189                            if (defined($essentialityData->{$temp[0]})) {
4190                                    if ($temp[1] eq "essential" && $essentialityData->{$temp[0]} eq "essential") {
4191                                            push(@{$predictionHash->{"Correct negative"}},$temp[0]);
4192                                    } elsif ($temp[1] eq "nonessential" && $essentialityData->{$temp[0]} eq "essential") {
4193                                            push(@{$predictionHash->{"False positive"}},$temp[0]);
4194                                    } elsif ($temp[1] eq "essential" && $essentialityData->{$temp[0]} eq "nonessential") {
4195                                            push(@{$predictionHash->{"False negative"}},$temp[0]);
4196                                    } elsif ($temp[1] eq "nonessential" && $essentialityData->{$temp[0]} eq "nonessential") {
4197                                            push(@{$predictionHash->{"Correct positive"}},$temp[0]);
4198                                    }
4199                            } else {
4200                                    push(@{$predictionHash->{$temp[1]}},$temp[0]);
4201                            }
4202                    }
4203                    foreach my $key (keys(%{$predictionHash})) {
4204                            my $string = $key;
4205                            $string = ucfirst($string);
4206                            push(@{$output},'<span title="'.join(", ",@{$predictionHash->{$key}}).'">'.$string.'</span>');
4207                    }
4208            }
4209    
4210            #Printing reactions
4211            if (defined($data->{$self->id()."REACTIONS"})) {
4212                    for (my $i=0; $i < @{$data->{$self->id()."REACTIONS"}};$i++) {
4213                            my $rxnData = $self->get_reaction_data($data->{$self->id()."REACTIONS"}->[$i]);
4214                            my $reactionString = $self->figmodel()->web()->create_reaction_link($data->{$self->id()."REACTIONS"}->[$i],join(" or ",@{$rxnData->{"ASSOCIATED PEG"}}),$self->id());
4215                            if (defined($rxnData->{PREDICTIONS})) {
4216                                    my $predictionHash;
4217                                    for (my $i=0; $i < @{$rxnData->{PREDICTIONS}};$i++) {
4218                                            my @temp = split(/:/,$rxnData->{PREDICTIONS}->[$i]);
4219                                            push(@{$predictionHash->{$temp[1]}},$temp[0]);
4220                                    }
4221                                    $reactionString .= "(";
4222                                    foreach my $key (keys(%{$predictionHash})) {
4223                                            if ($key eq "Essential =>") {
4224                                                    $reactionString .= '<span title="Essential in '.join(",",@{$predictionHash->{$key}}).'">E=></span>,';
4225                                            } elsif ($key eq "Essential <=") {
4226                                                    $reactionString .= '<span title="Essential in '.join(",",@{$predictionHash->{$key}}).'">E<=</span>,';
4227                                            } elsif ($key eq "Active =>") {
4228                                                    $reactionString .= '<span title="Active in '.join(",",@{$predictionHash->{$key}}).'">A=></span>,';
4229                                            } elsif ($key eq "Active <=") {
4230                                                    $reactionString .= '<span title="Active in '.join(",",@{$predictionHash->{$key}}).'">A<=</span>,';
4231                                            } elsif ($key eq "Active <=>") {
4232                                                    $reactionString .= '<span title="Active in '.join(",",@{$predictionHash->{$key}}).'">A</span>,';
4233                                            } elsif ($key eq "Inactive") {
4234                                                    $reactionString .= '<span title="Inactive in '.join(",",@{$predictionHash->{$key}}).'">I</span>,';
4235                                            } elsif ($key eq "Dead") {
4236                                                    $reactionString .= '<span title="Dead">D</span>,';
4237                                            }
4238                                    }
4239                                    $reactionString =~ s/,$/)/;
4240                            }
4241                            push(@{$output},$reactionString);
4242                    }
4243            }
4244    
4245            #Returning output
4246            return join("<br>",@{$output});
4247    }
4248    
4249    =head3 remove_obsolete_reactions
4250    Definition:
4251            void FIGMODELmodel->remove_obsolete_reactions();
4252    Description:
4253    =cut
4254    sub remove_obsolete_reactions {
4255            my ($self) = @_;
4256    
4257            (my $dummy,my $translation) = $self->figmodel()->put_two_column_array_in_hash($self->figmodel()->database()->load_multiple_column_file($self->figmodel()->config("Translation directory")->[0]."ObsoleteRxnIDs.txt","\t"));
4258            my $rxnTbl = $self->reaction_table();
4259            if (defined($rxnTbl)) {
4260                    for (my $i=0; $i < $rxnTbl->size(); $i++) {
4261                            my $row = $rxnTbl->get_row($i);
4262                            if (defined($translation->{$row->{LOAD}->[0]}) || defined($translation->{$row->{LOAD}->[0]."r"})) {
4263                                    my $direction = $row->{DIRECTION}->[0];
4264                                    my $newRxn;
4265                                    if (defined($translation->{$row->{LOAD}->[0]."r"})) {
4266                                            $newRxn = $translation->{$row->{LOAD}->[0]."r"};
4267                                            if ($direction eq "<=") {
4268                                                    $direction = "=>";
4269                                            } elsif ($direction eq "=>") {
4270                                                    $direction = "<=";
4271                                            }
4272                                    } else {
4273                                            $newRxn = $translation->{$row->{LOAD}->[0]};
4274                                    }
4275                                    #Checking if the new reaction is already in the model
4276                                    my $newRow = $rxnTbl->get_row_by_key($newRxn,"LOAD");
4277                                    if (defined($newRow)) {
4278                                            #Handling direction
4279                                            if ($newRow->{DIRECTION}->[0] ne $direction) {
4280                                                    $newRow->{DIRECTION}->[0] = "<=>";
4281                                            }
4282                                            push(@{$row->{"ASSOCIATED PEG"}},@{$rxnTbl->get_row($i)->{"ASSOCIATED PEG"}});
4283                                    } else {
4284                                            $rxnTbl->get_row($i)->{LOAD}->[0] = $newRxn;
4285                                            $rxnTbl->get_row($i)->{DIRECTION}->[0] = $direction;
4286                                    }
4287                            }
4288                    }
4289                    $rxnTbl->save();
4290            }
4291    }
4292    
4293  1;  1;

Legend:
Removed from v.1.2  
changed lines
  Added in v.1.14

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3