[Bio] / Sprout / GenomeStats.pl Repository:
ViewVC logotype

Diff of /Sprout/GenomeStats.pl

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

revision 1.4, Sun Jun 18 07:12:58 2006 UTC revision 1.30, Thu Dec 6 14:58:03 2007 UTC
# Line 6  Line 6 
6  the genomes in each of the genome groups. Genomes that are new to this version  the genomes in each of the genome groups. Genomes that are new to this version
7  of the Sprout will be specially marked. In order for this to work, both the  of the Sprout will be specially marked. In order for this to work, both the
8  current and previous Sprout databases must be available on this machine.  current and previous Sprout databases must be available on this machine.
 This is one positional parameter: the name of a directory in which to place  
 the include files.  
9    
10  The currently-supported command-line options are as follows.  The currently-supported command-line options are as follows.
11    
# Line 61  Line 59 
59    
60  Style to use for small-text markers (e.g. NEW!)  Style to use for small-text markers (e.g. NEW!)
61    
62    =item numStyle
63    
64    Style to use for numeric cells.
65    
66    =item counterStyle
67    
68    Style to use for counter cells.
69    
70  =item linkCGI  =item linkCGI
71    
72  Path to the CGI script for displaying detailed statistics.  Path to the CGI script for displaying detailed statistics.
73    
74    =item noNewCheck
75    
76    If specified, the check for new genomes in the group is suppressed. This
77    may need to be done if there's been a change in the database definition. Note
78    that all this really does is keep the B<NEW!> symbol from showing. It does
79    not affect which genomes show up in the table.
80    
81  =back  =back
82    
83  =cut  =cut
# Line 80  Line 93 
93  use SFXlate;  use SFXlate;
94  use CGI qw(:standard);  use CGI qw(:standard);
95  use FIG;  use FIG;
 no warnings 'once'; # only when coding  
96    
97  # Get the command-line options and parameters.  # Get the command-line options and parameters.
98  my ($options, @parameters) = StandardSetup([qw(Sprout ERDB) ],  my ($options, @parameters) = StandardSetup([qw(Sprout ERDB) ],
99                                             {                                             {
100                                              strict => [0, 'keep related groups separate'],                                              strict => [0, 'keep related groups separate'],
101                                              oddStyle => ['odd', 'style for odd rows'],                                              oddStyle => ['odd', 'style for odd rows'],
102                                                trace => [2, 'tracing level'],
103                                              evenStyle => ['even', 'style for even rows'],                                              evenStyle => ['even', 'style for even rows'],
104                                              tableStyle => ['genomestats', 'style for whole table'],                                              tableStyle => ['genomestats', 'style for whole table'],
105                                              markerStyle => ['tinytext', 'style for markers'],                                              markerStyle => ['tinytext', 'style for markers'],
106                                                numStyle => ['numcell', 'style for cells with numeric values'],
107                                                counterStyle => ['countercell', 'style for cells with counter values'],
108                                              linkCGI => ['../FIG/genome_statistics.cgi',                                              linkCGI => ['../FIG/genome_statistics.cgi',
109                                                          'path to CGI script for detailed statistics'],                                                          'path to CGI script for detailed statistics'],
110                                                noNewCheck => [0, 'if specified, skips the check for new genomes'],
111                                                targetDir => ["$FIG_Config::nmpdr_base/next/html/includes",
112                                                              'target directory'],
113                                             },                                             },
114                                             "<targetDir>",                                             "",
115                                             @ARGV);                                             @ARGV);
116    # The return type (error/no error) goes in here.
117    my $rtype;
118    eval {
119        # This table controls the special attribute columns. For each we need to know the attribute name and the
120        # column title. If any genomes in a group have a value for one of the special columns, that column is
121        # displayed along with the attribute values.
122        my %specialCols = (Serotype => 'Serotype_code',
123                           Phenotype => 'Phenotype');
124  # Verify the directory name.  # Verify the directory name.
125  my $targetDir = $parameters[0];      my $targetDir = $options->{targetDir};
126  if (! $targetDir) {  if (! $targetDir) {
127      Confess("No target directory specified.");      Confess("No target directory specified.");
128  } elsif (! -d $targetDir) {  } elsif (! -d $targetDir) {
129      Confess("Target directory $targetDir not found.");      Confess("Target directory $targetDir not found.");
130  } else {  } else {
     # *Get the old Sprout.  
     my $oldSprout = SFXlate->new_sprout_only($FIG_Config::oldSproutDB);  
     # Extract the genome group data from the old Sprout.  
     my %oldGroupHash = $oldSprout->GetGroups();  
     if (! $options->{strict}) {  
         %oldGroupHash = Fix(%oldGroupHash);  
     }  
131      # Get the new Sprout.      # Get the new Sprout.
132      my $sprout = SFXlate->new_sprout_only();      my $sprout = SFXlate->new_sprout_only();
133      my %newGroupHash = $sprout->GetGroups();      my %newGroupHash = $sprout->GetGroups();
134            # Extract the genome group data from the new Sprout.
135      if (! $options->{strict}) {      if (! $options->{strict}) {
136          %newGroupHash = Fix(%newGroupHash);              %newGroupHash = $sprout->Fix(%newGroupHash);
137      }      }
138            # This hash will be used to determine which genomes are new.
139            my %oldGroupHash = ();
140            if ($options->{noNewCheck}) {
141                # Here we can't look at the old Sprout. Set up the hash
142                # so it looks like the old Sprout's data is the same as ours.
143                %oldGroupHash = map { $_ => $newGroupHash{$_} } keys %newGroupHash;
144            } else {
145                # Get the old Sprout.
146                my $oldSprout = SFXlate->old_sprout_only();
147                # Extract the genome group data from the old Sprout.
148                %oldGroupHash = $oldSprout->GetGroups();
149                if (! $options->{strict}) {
150                    %oldGroupHash = $oldSprout->Fix(%oldGroupHash);
151                }
152            }
153            # Get a FIG object for computing attributes.
154            my $fig = FIG->new();
155            # Get the super-group list.
156            my @superGroups = sort keys %newGroupHash;
157            # Set up some useful stuff for the four count columns.
158            my %linkParms = ( s0 => "nothypo_sub", n0 => "nothypo_nosub",
159                              s1 => "hypo_sub", n1 => "hypo_nosub" );
160            my @columnTypes = ('s0', 'n0', 's1', 'n1');
161            # Get the styles.
162            my ($tableStyle, $markerStyle, @rowStyle) = ($options->{tableStyle}, $options->{markerStyle},
163                                                         $options->{evenStyle}, $options->{oddStyle});
164            my ($numStyle, $counterStyle) = ($options->{numStyle}, $options->{counterStyle});
165            # Prepare a hash for the summary counters. These will be used on the organism summary page.
166            my %summaries = ();
167      # Loop through the groups.      # Loop through the groups.
168      for my $groupID (keys %newGroupHash) {          for my $groupID (@superGroups) {
169          Trace("Processing group $groupID.") if T(2);          Trace("Processing group $groupID.") if T(2);
170                # Create a hash for summarizing the counters.
171                my %groupTotals = ( genomes => 0, pegs => 0, RNAs => 0,
172                                   map { $_ => } @columnTypes, features => 0 );
173          # Get the genomes from the new hash.          # Get the genomes from the new hash.
174          my @newGenomes = @{$newGroupHash{$groupID}};          my @newGenomes = @{$newGroupHash{$groupID}};
175          # Create a hash for finding if a genome is in the old group. If the entire group is          # Create a hash for finding if a genome is in the old group. If the entire group is
# Line 127  Line 179 
179              %oldGenomes = map { $_ => 1 } @{$oldGroupHash{$groupID}};              %oldGenomes = map { $_ => 1 } @{$oldGroupHash{$groupID}};
180          }          }
181          # Create the output file.          # Create the output file.
182          Open(\*GROUPFILE, ">$targetDir/$groupID.inc");              my $outFileName = "stats-" . lc($groupID) . ".inc";
183          # Get the styles.              Open(\*GROUPFILE, ">$targetDir/$outFileName");
184          my ($tableStyle, $markerStyle, @rowStyle) = ($options->{tableStyle}, $options->{markerStyle},              # Get the special columns. We'll stuff them in a hash keyed by column name. Each column name will contain
185                                                       $options->{evenStyle}, $options->{oddStyle});              # a sub-hash that translates each genome ID to its applicable attribute value (if any).
186          # Start the table.              my %specialData = ();
187          print GROUPFILE "<table class=\"$tableStyle\">\n";              for my $specialColumn (keys %specialCols) {
188          # Create the header row.                  # Get the attribute mapping.
189          print GROUPFILE Tr( { class => 'odd' }, th("Strain annotated in NMPDR",                  my %specialDataList = map { $_->[0] => $_->[2] } $fig->get_attributes(\@newGenomes, $specialCols{$specialColumn});
190                                                   "Genome size, bp",                  # We only proceed if some attributes were found. As a result, the keys in %specialData will only be keys
191                    # for columns that exist in the output.
192                    if (scalar(keys %specialDataList)) {
193                        $specialData{$specialColumn} = \%specialDataList;
194                    }
195                }
196                # Set up the column names.
197                my @columnNames = "Strain annotated in NMPDR";
198                push @columnNames, sort keys %specialData;
199                push @columnNames,  "Genome size, bp",
200                                                   "Protein Encoding Genes (PEGs)",                                                   "Protein Encoding Genes (PEGs)",
201                                                   "Named genes in subsystems",            # s0                                                   "Named genes in subsystems",            # s0
202                                                   "Named genes not in subsystems",        # n0                                                   "Named genes not in subsystems",        # n0
203                                                   "Hypothetical genes in subsystems",     # s1                                                   "Hypothetical genes in subsystems",     # s1
204                                                   "Hypothetical genes not in subsystems", # n1                                                   "Hypothetical genes not in subsystems", # n1
205                                                   "RNAs")) . "\n";                                  "Subsystems",
206          # Set up some useful stuff for the four count columns.                                  "RNAs";
207          my %linkParms = ( s0 => "nohypo_sub", n0 => "nohypo_nosub",              # Start the table.
208                            s1 => "hypo_sub", n1 => "hypo_nosub" );              print GROUPFILE "<table class=\"$tableStyle\">\n";
209          my @columnTypes = ('s0', 'n0', 's1', 'n1');              # Create the header row.
210                print GROUPFILE Tr( { class => 'odd' }, th(\@columnNames)) . "\n";
211          # The data rows will be built next. We'll be putting them into a hash keyed by          # The data rows will be built next. We'll be putting them into a hash keyed by
212          # organism name. The hash enables us to spit them out sorted by name.          # organism name. The hash enables us to spit them out sorted by name.
213          my %rows = ();          my %rows = ();
214                # This variable is used to hold the counts.
215                my $num;
216          # Loop through the genomes in the new group.          # Loop through the genomes in the new group.
217          for my $genomeID (@newGenomes) {          for my $genomeID (@newGenomes) {
218                    # Count this genome.
219                    $groupTotals{genomes}++;
220              # Check to see if this genome is new.              # Check to see if this genome is new.
221              my $new = (! exists $oldGenomes{$genomeID} ? "new " : "");              my $new = (! exists $oldGenomes{$genomeID} ? "new " : "");
222              Trace("Processing ${new}genome $genomeID for $groupID.") if T(3);              Trace("Processing ${new}genome $genomeID for $groupID.") if T(3);
223              # Get the strain name.              # Get the strain name.
224              my $genomeName = $sprout->GenusSpecies($genomeID);              my $genomeName = $sprout->GenusSpecies($genomeID);
225                    # Apply a link.
226                    my $genomeText = CGI::a({ href => "../FIG/genome_statistics.cgi?genome=$genomeID;SPROUT=1" }, $genomeName);
227              # If this is a new strain, build the HTML for the NEW! mark.              # If this is a new strain, build the HTML for the NEW! mark.
228              if ($new) {              if ($new) {
229                  $new = " <span class=\"$markerStyle\">NEW!</span>";                  $new = " <span class=\"$markerStyle\">NEW!</span>";
230              }              }
231              # Get the genome length.              # Get the genome length.
232              my $genomeLen = $sprout->GenomeLength($genomeID);                  $num = $sprout->GenomeLength($genomeID);
233                    my $genomeLen = Tracer::CommaFormat($num);
234              # Get the number of PEGs.              # Get the number of PEGs.
235              my $pegCount = $sprout->FeatureCount($genomeID, 'peg');                  $num = $sprout->FeatureCount($genomeID, 'peg');
236                    my $pegCount = Tracer::CommaFormat($num);
237                    $groupTotals{pegs} += $num;
238              # Get the number of RNAs.              # Get the number of RNAs.
239              my $rnaCount = $sprout->FeatureCount($genomeID, 'rna');                  $num = $sprout->FeatureCount($genomeID, 'rna');
240                    my $rnaCount = Tracer::CommaFormat($num);
241                    $groupTotals{RNAs} += $num;
242                    # If there are no RNAs, we say we don't know the number, since we know there
243                    # must be RNAs somewhere.
244                    if (! $rnaCount) {
245                        $rnaCount = "n/d";
246                    }
247              # Now we have four categories of features to work with, for each              # Now we have four categories of features to work with, for each
248              # combination of named or hypothetical vs. in-subsystem or              # combination of named or hypothetical vs. in-subsystem or
249              # not-in-subsystem. First, we get all of the feature assignments for              # not-in-subsystem. First, we get all of the feature assignments for
250              # the genome.              # the genome.
251              my $assignHash = $sprout->GenomeAssignments($genomeID);              my $assignHash = $sprout->GenomeAssignments($genomeID);
252              # Next, we get all of the features in the genome that belong to a              # Next, we get all of the features in the genome that belong to a
253              # subsystem. This involves a query via the subsystem spreadsheet.                  # subsystem.
254              my %ssHash = map { $_ => 1 } $sprout->GetFlat(['IsGenomeOf', 'ContainsFeature'],                  my %ssHash = $sprout->GenomeSubsystemData($genomeID);
                                                     "IsGenomeOf(from-link) = ?",  
                                                     [$genomeID], 'ContainsFeature(to-link)');  
255              # Create a hash to track the four categories. "s" or "n" indicates              # Create a hash to track the four categories. "s" or "n" indicates
256              # in or out of a subsystem. "1" or "0" indicates hypothetical or              # in or out of a subsystem. "1" or "0" indicates hypothetical or
257              # real.              # real.
# Line 190  Line 266 
266                  $counters{$ss} += 1;                  $counters{$ss} += 1;
267                  $totalFeatures++;                  $totalFeatures++;
268              }              }
269              Trace("$totalFeatures total feature found for $genomeID.") if T(3);                  Trace("$totalFeatures total features found for $genomeID.") if T(3);
270                    for my $counterKey (@columnTypes) {
271                        $groupTotals{$counterKey} += $counters{$counterKey};
272                    }
273                    $groupTotals{features} += $totalFeatures;
274              # We have all our data. Next we need to compute the percentages and the links.              # We have all our data. Next we need to compute the percentages and the links.
275              # First, the link stuff.              # First, the link stuff.
276              my $linkPrefix = "$options->{linkCGI}?user=\&genome=$genomeID&SPROUT=1&request=";              my $linkPrefix = "$options->{linkCGI}?user=\&genome=$genomeID&SPROUT=1&request=";
# Line 198  Line 278 
278              for my $type (keys %linkParms) {              for my $type (keys %linkParms) {
279                  $counters{$type} = a( { href => "$linkPrefix$linkParms{$type}" },                  $counters{$type} = a( { href => "$linkPrefix$linkParms{$type}" },
280                                       sprintf("%d(%.1f%%)", $counters{$type},                                       sprintf("%d(%.1f%%)", $counters{$type},
281                                               $counters{$type} * 100 / $totalFeatures));                                                   Tracer::Percent($counters{$type}, $totalFeatures)));
282              }              }
283              # Create the row text.                  my @counterValues = map { $counters{$_} } @columnTypes;
284              my $rowHtml = td( "$genomeName$new", $genomeLen, $pegCount,                  # The last link is a button to look at the subsystem summaries.
285                                map { $counters{$_} } @columnTypes,                  my $ssCount = $sprout->GetCount(['ParticipatesIn'], 'ParticipatesIn(from-link) = ?',
286                                $rnaCount );                                                     [$genomeID]);
287                    my $ssLink = "$options->{linkCGI}?user=\&genome=$genomeID&SPROUT=1&show_subsystems=1";
288                    my $ssCol = "<a href=\"$ssLink\">$ssCount</a>";
289                    # Start creating the table cells.
290                    my $rowHtml = td("$genomeText$new");
291                    # Add any special columns.
292                    for my $specialCol (keys %specialData) {
293                        # Here we get the attribute value. If there is none, we leave the column blank.
294                        my $attribute = $specialData{$specialCol}->{$genomeID} || "&nbsp;";
295                        $rowHtml .= td($attribute);
296                    }
297                    # Now add the data columns.
298                    $rowHtml .= join("",
299                                    td({ class => $numStyle }, $genomeLen),
300                                    td({ class => $numStyle }, $pegCount),
301                                    td({ class => $counterStyle }, \@counterValues),
302                                    td({ class => $numStyle }, $ssCol),
303                                    td({ class => $numStyle }, $rnaCount),
304                                    );
305              # Put it in the row hash.              # Put it in the row hash.
306              $rows{$genomeName} = $rowHtml;              $rows{$genomeName} = $rowHtml;
307          }          }
# Line 222  Line 320 
320              # Count the row.              # Count the row.
321              $rowCount++;              $rowCount++;
322          }          }
323          # All done, close the file.              # All done, terminate the table and close the file.
324                print GROUPFILE "</table>\n";
325          close GROUPFILE;          close GROUPFILE;
326          Trace("$rowCount genomes processed.") if T(2);          Trace("$rowCount genomes processed.") if T(2);
327                # Now save the group totals.
328                $summaries{$groupID} = \%groupTotals;
329      }      }
330            # Now produce the summary table.
331            my $sumFileName = "stats-groups.inc";
332            Open(\*SUMFILE, ">$targetDir/$sumFileName");
333            # Start the table.
334            print SUMFILE "<table class=\"$tableStyle\">\n";
335            # Create the header row.
336            print SUMFILE Tr( { class => 'odd' }, th(["Group name",
337                                                     "Genomes",
338                                                     "Protein Encoding Genes (PEGs)",
339                                                     "Named genes in subsystems",            # s0
340                                                     "Named genes not in subsystems",        # n0
341                                                     "Hypothetical genes in subsystems",     # s1
342                                                     "Hypothetical genes not in subsystems", # n1
343                                                     "RNAs",
344                                                       ])) . "\n";
345            # Set up a flag for the odd-even styling.
346            my $rowFlag = 0;
347            # Put in the data rows.
348            for my $groupName (sort keys %summaries) {
349                my $group = $summaries{$groupName};
350                # Compute the link for the current group.
351                my $groupLink = a({ href => $sprout->GroupPageName($groupName) }, $groupName);
352                # Create the table row.
353                my $rowHtml = join("",
354                                   td($groupLink),
355                                   td({ class => $numStyle }, Tracer::CommaFormat($group->{genomes})),
356                                   td({ class => $numStyle }, Tracer::CommaFormat($group->{pegs})),
357                                   td({ class => $counterStyle }, [ map { Tracer::CommaFormat($group->{$_}) } @columnTypes ]),
358                                   td({ class => $numStyle }, Tracer::CommaFormat($group->{RNAs})),
359                                  );
360                print SUMFILE Tr( { class => $rowStyle[$rowFlag] }, $rowHtml ) . "\n";
361                # Flip the row style.
362                $rowFlag = 1 - $rowFlag;
363            }
364            # Terminate the table and close the file.
365            print SUMFILE "</table>\n";
366            close SUMFILE;
367      # We're all done.      # We're all done.
368      Trace("Processing complete.") if T(2);      Trace("Processing complete.") if T(2);
369  }  }
370    };
371  =head3 Fix  if ($@) {
372        Trace("Stats failed with error: $@") if T(0);
373  C<< my %fixedHash = Fix(%groupHash); >>      $rtype = "error";
374    } else {
375  Prepare a genome group hash for processing. Groups with the same primary name will be combined.      Trace("Stats complete.") if T(2);
376  The primary name is the first capitalized word in the group name.      $rtype = "no error";
   
 =over 4  
   
 =item groupHash  
   
 Hash to be fixed up.  
   
 =item RETURN  
   
 Returns a fixed-up version of the hash.  
   
 =back  
   
 =cut  
   
 sub Fix {  
     # Get the parameters.  
     my (%groupHash) = @_;  
     # Create the result hash.  
     my %retVal = ();  
     # Copy over the genomes.  
     for my $groupID (keys %groupHash) {  
         # Make a safety copy of the group ID.  
         my $realGroupID = $groupID;  
         # Yank the primary name.  
         if ($groupID =~ /([A-Z]\w+)/) {  
             $realGroupID = $1;  
377          }          }
378          # Append this group's genomes into the result hash.  if ($options->{phone}) {
379          Tracer::AddToListMap(\%retVal, $realGroupID, $groupHash{$groupID});      my $msgID = Tracer::SendSMS($options->{phone}, "GenomeStats terminated with $rtype.");
380        if ($msgID) {
381            Trace("Phone message sent with ID $msgID.") if T(2);
382        } else {
383            Trace("Phone message not sent.") if T(2);
384      }      }
     # Return the result hash.  
     return %retVal;  
385  }  }
386    
   
387  1;  1;
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.30
MCS Webmaster
 ViewVC Help
Powered by ViewVC 1.0.3