[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.15, Fri Jun 23 23:12:43 2006 UTC revision 1.28, Tue Apr 10 06:01:56 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 73  Line 71 
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 95  Line 100 
100                                             {                                             {
101                                              strict => [0, 'keep related groups separate'],                                              strict => [0, 'keep related groups separate'],
102                                              oddStyle => ['odd', 'style for odd rows'],                                              oddStyle => ['odd', 'style for odd rows'],
103                                                trace => [2, 'tracing level'],
104                                              evenStyle => ['even', 'style for even rows'],                                              evenStyle => ['even', 'style for even rows'],
105                                              tableStyle => ['genomestats', 'style for whole table'],                                              tableStyle => ['genomestats', 'style for whole table'],
106                                              markerStyle => ['tinytext', 'style for markers'],                                              markerStyle => ['tinytext', 'style for markers'],
# Line 102  Line 108 
108                                              counterStyle => ['countercell', 'style for cells with counter values'],                                              counterStyle => ['countercell', 'style for cells with counter values'],
109                                              linkCGI => ['../FIG/genome_statistics.cgi',                                              linkCGI => ['../FIG/genome_statistics.cgi',
110                                                          'path to CGI script for detailed statistics'],                                                          'path to CGI script for detailed statistics'],
111                                                groupFile => ["$FIG_Config::sproutData/groups.tbl",
112                                                              'location of the NMPDR group description file'],
113                                                noNewCheck => [0, 'if specified, skips the check for new genomes'],
114                                                targetDir => ["$FIG_Config::nmpdr_base/next/html/includes",
115                                                              'target directory'],
116                                             },                                             },
117                                             "<targetDir>",                                             "",
118                                             @ARGV);                                             @ARGV);
119  # Verify the directory name.  # Verify the directory name.
120  my $targetDir = $parameters[0];  my $targetDir = $options->{targetDir};
121  if (! $targetDir) {  if (! $targetDir) {
122      Confess("No target directory specified.");      Confess("No target directory specified.");
123  } elsif (! -d $targetDir) {  } elsif (! -d $targetDir) {
124      Confess("Target directory $targetDir not found.");      Confess("Target directory $targetDir not found.");
125  } 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);  
     }  
126      # Get the new Sprout.      # Get the new Sprout.
127      my $sprout = SFXlate->new_sprout_only();      my $sprout = SFXlate->new_sprout_only();
128      my %newGroupHash = $sprout->GetGroups();      my %newGroupHash = $sprout->GetGroups();
129        # Extract the genome group data from the new Sprout.
130        if (! $options->{strict}) {
131            %newGroupHash = Sprout::Fix(%newGroupHash);
132        }
133        # This hash will be used to determine which genomes are new.
134        my %oldGroupHash = ();
135        if ($options->{noNewCheck}) {
136            # Here we can't look at the old Sprout. Set up the hash
137            # so it looks like the old Sprout's data is the same as ours.
138            %oldGroupHash = map { $_ => $newGroupHash{$_} } keys %newGroupHash;
139        } else {
140            # Get the old Sprout.
141            my $oldSprout = SFXlate->old_sprout_only();
142            # Extract the genome group data from the old Sprout.
143            %oldGroupHash = $oldSprout->GetGroups();
144      if (! $options->{strict}) {      if (! $options->{strict}) {
145          %newGroupHash = Fix(%newGroupHash);              %oldGroupHash = Sprout::Fix(%oldGroupHash);
146      }      }
147        }
148        # Get a FIG object for computing attributes.
149        my $fig = FIG->new();
150        # Read the group file.
151        my %groupData = Sprout::ReadGroupFile($options->{groupFile});
152        # Set up some useful stuff for the four count columns.
153        my %linkParms = ( s0 => "nothypo_sub", n0 => "nothypo_nosub",
154                          s1 => "hypo_sub", n1 => "hypo_nosub" );
155        my @columnTypes = ('s0', 'n0', 's1', 'n1');
156        # Get the styles.
157        my ($tableStyle, $markerStyle, @rowStyle) = ($options->{tableStyle}, $options->{markerStyle},
158                                                     $options->{evenStyle}, $options->{oddStyle});
159        my ($numStyle, $counterStyle) = ($options->{numStyle}, $options->{counterStyle});
160        # Prepare a hash for the summary counters. These will be used on the organism summary page.
161        my %summaries = ();
162      # Loop through the groups.      # Loop through the groups.
163      for my $groupID (keys %newGroupHash) {      for my $groupID (keys %newGroupHash) {
164          Trace("Processing group $groupID.") if T(2);          Trace("Processing group $groupID.") if T(2);
165            # Create a hash for summarizing the counters.
166            my %groupTotals = ( genomes => 0, pegs => 0, RNAs => 0,
167                               map { $_ => } @columnTypes, features => 0 );
168          # Get the genomes from the new hash.          # Get the genomes from the new hash.
169          my @newGenomes = @{$newGroupHash{$groupID}};          my @newGenomes = @{$newGroupHash{$groupID}};
170          # 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 139  Line 176 
176          # Create the output file.          # Create the output file.
177          my $outFileName = "stats-" . lc($groupID) . ".inc";          my $outFileName = "stats-" . lc($groupID) . ".inc";
178          Open(\*GROUPFILE, ">$targetDir/$outFileName");          Open(\*GROUPFILE, ">$targetDir/$outFileName");
179          # Get the styles.          # Get the serotypes.
180          my ($tableStyle, $markerStyle, @rowStyle) = ($options->{tableStyle}, $options->{markerStyle},          my %serotypes = map { $_->[0] => $_->[2] } $fig->get_attributes(\@newGenomes, "Serotype_code");
181                                                       $options->{evenStyle}, $options->{oddStyle});          my $hasSeroData = (scalar(keys %serotypes) > 0);
182          my ($numStyle, $counterStyle) = ($options->{numStyle}, $options->{counterStyle});          # If we have serotypes, we add an extra column.
183          # Start the table.          my @columnNames = "Strain annotated in NMPDR";
184          print GROUPFILE "<table class=\"$tableStyle\">\n";          if ($hasSeroData) {
185          # Create the header row.              push @columnNames, "Serotype";
186          print GROUPFILE Tr( { class => 'odd' }, th(["Strain annotated in NMPDR",          }
187                                                   "Genome size, bp",          push @columnNames,  "Genome size, bp",
188                                                   "Protein Encoding Genes (PEGs)",                                                   "Protein Encoding Genes (PEGs)",
189                                                   "Named genes in subsystems",            # s0                                                   "Named genes in subsystems",            # s0
190                                                   "Named genes not in subsystems",        # n0                                                   "Named genes not in subsystems",        # n0
191                                                   "Hypothetical genes in subsystems",     # s1                                                   "Hypothetical genes in subsystems",     # s1
192                                                   "Hypothetical genes not in subsystems", # n1                                                   "Hypothetical genes not in subsystems", # n1
193                                                   "RNAs",                              "Subsystems",
194                                                     ])) . "\n";                              "RNAs";
195          # Set up some useful stuff for the four count columns.          # Start the table.
196          my %linkParms = ( s0 => "nothypo_sub", n0 => "nothypo_nosub",          print GROUPFILE "<table class=\"$tableStyle\">\n";
197                            s1 => "hypo_sub", n1 => "hypo_nosub" );          # Create the header row.
198          my @columnTypes = ('s0', 'n0', 's1', 'n1');          print GROUPFILE Tr( { class => 'odd' }, th(\@columnNames)) . "\n";
199          # 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
200          # 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.
201          my %rows = ();          my %rows = ();
202            # This variable is used to hold the counts.
203            my $num;
204          # Loop through the genomes in the new group.          # Loop through the genomes in the new group.
205          for my $genomeID (@newGenomes) {          for my $genomeID (@newGenomes) {
206                # Count this genome.
207                $groupTotals{genomes}++;
208              # Check to see if this genome is new.              # Check to see if this genome is new.
209              my $new = (! exists $oldGenomes{$genomeID} ? "new " : "");              my $new = (! exists $oldGenomes{$genomeID} ? "new " : "");
210              Trace("Processing ${new}genome $genomeID for $groupID.") if T(3);              Trace("Processing ${new}genome $genomeID for $groupID.") if T(3);
211              # Get the strain name.              # Get the strain name.
212              my $genomeName = $sprout->GenusSpecies($genomeID);              my $genomeName = $sprout->GenusSpecies($genomeID);
213                # Apply a link.
214                my $genomeText = CGI::a({ href => "../FIG/genome_statistics.cgi?genome=$genomeID;SPROUT=1" }, $genomeName);
215              # 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.
216              if ($new) {              if ($new) {
217                  $new = " <span class=\"$markerStyle\">NEW!</span>";                  $new = " <span class=\"$markerStyle\">NEW!</span>";
218              }              }
219              # Get the genome length.              # Get the genome length.
220              my $genomeLen = Tracer::CommaFormat($sprout->GenomeLength($genomeID));              $num = $sprout->GenomeLength($genomeID);
221                my $genomeLen = Tracer::CommaFormat($num);
222              # Get the number of PEGs.              # Get the number of PEGs.
223              my $pegCount = Tracer::CommaFormat($sprout->FeatureCount($genomeID, 'peg'));              $num = $sprout->FeatureCount($genomeID, 'peg');
224                my $pegCount = Tracer::CommaFormat($num);
225                $groupTotals{pegs} += $num;
226              # Get the number of RNAs.              # Get the number of RNAs.
227              my $rnaCount = Tracer::CommaFormat($sprout->FeatureCount($genomeID, 'rna'));              $num = $sprout->FeatureCount($genomeID, 'rna');
228                my $rnaCount = Tracer::CommaFormat($num);
229                $groupTotals{RNAs} += $num;
230                # If there are no RNAs, we say we don't know the number, since we know there
231                # must be RNAs somewhere.
232                if (! $rnaCount) {
233                    $rnaCount = "n/d";
234                }
235              # Now we have four categories of features to work with, for each              # Now we have four categories of features to work with, for each
236              # combination of named or hypothetical vs. in-subsystem or              # combination of named or hypothetical vs. in-subsystem or
237              # not-in-subsystem. First, we get all of the feature assignments for              # not-in-subsystem. First, we get all of the feature assignments for
238              # the genome.              # the genome.
239              my $assignHash = $sprout->GenomeAssignments($genomeID);              my $assignHash = $sprout->GenomeAssignments($genomeID);
240              # 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
241              # subsystem. This involves a query via the subsystem spreadsheet.              # subsystem.
242              my %ssHash = map { $_ => 1 } $sprout->GetFlat(['IsGenomeOf', 'ContainsFeature'],              my %ssHash = $sprout->GenomeSubsystemData($genomeID);
                                                     "IsGenomeOf(from-link) = ?",  
                                                     [$genomeID], 'ContainsFeature(to-link)');  
243              # Create a hash to track the four categories. "s" or "n" indicates              # Create a hash to track the four categories. "s" or "n" indicates
244              # in or out of a subsystem. "1" or "0" indicates hypothetical or              # in or out of a subsystem. "1" or "0" indicates hypothetical or
245              # real.              # real.
# Line 204  Line 255 
255                  $totalFeatures++;                  $totalFeatures++;
256              }              }
257              Trace("$totalFeatures total features found for $genomeID.") if T(3);              Trace("$totalFeatures total features found for $genomeID.") if T(3);
258                for my $counterKey (@columnTypes) {
259                    $groupTotals{$counterKey} += $counters{$counterKey};
260                }
261                $groupTotals{features} += $totalFeatures;
262              # 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.
263              # First, the link stuff.              # First, the link stuff.
264              my $linkPrefix = "$options->{linkCGI}?user=\&genome=$genomeID&SPROUT=1&request=";              my $linkPrefix = "$options->{linkCGI}?user=\&genome=$genomeID&SPROUT=1&request=";
# Line 214  Line 269 
269                                               Tracer::Percent($counters{$type}, $totalFeatures)));                                               Tracer::Percent($counters{$type}, $totalFeatures)));
270              }              }
271              my @counterValues = map { $counters{$_} } @columnTypes;              my @counterValues = map { $counters{$_} } @columnTypes;
272              # Create the row text. Note that we use the distributive capability of the TD              # The last link is a button to look at the subsystem summaries.
273              # function to apply the same style to each one.              my $ssCount = $sprout->GetCount(['ParticipatesIn'], 'ParticipatesIn(from-link) = ?',
274              my $rowHtml = join("",                                                 [$genomeID]);
275                                 td("$genomeName$new"),              my $ssLink = "$options->{linkCGI}?user=\&genome=$genomeID&SPROUT=1&show_subsystems=1";
276                my $ssCol = "<a href=\"$ssLink\">$ssCount</a>";
277                # Start creating the table cells.
278                my $rowHtml = td("$genomeText$new");
279                # Check for a serotype.
280                if ($hasSeroData) {
281                    my $seroType = $serotypes{$genomeID} || "&nbsp;";
282                    $rowHtml .= td($seroType);
283                }
284                # Now add the data columns.
285                $rowHtml .= join("",
286                                 td({ class => $numStyle }, $genomeLen),                                 td({ class => $numStyle }, $genomeLen),
287                                 td({ class => $numStyle }, $pegCount),                                 td({ class => $numStyle }, $pegCount),
288                                 td({ class => $counterStyle }, \@counterValues),                                 td({ class => $counterStyle }, \@counterValues),
289                                td({ class => $numStyle }, $ssCol),
290                                 td({ class => $numStyle }, $rnaCount),                                 td({ class => $numStyle }, $rnaCount),
291                                );                                );
292              # Put it in the row hash.              # Put it in the row hash.
# Line 245  Line 311 
311          print GROUPFILE "</table>\n";          print GROUPFILE "</table>\n";
312          close GROUPFILE;          close GROUPFILE;
313          Trace("$rowCount genomes processed.") if T(2);          Trace("$rowCount genomes processed.") if T(2);
314            # Now save the group totals.
315            $summaries{$groupID} = \%groupTotals;
316      }      }
317        # Now produce the summary table.
318        my $sumFileName = "stats-groups.inc";
319        Open(\*SUMFILE, ">$targetDir/$sumFileName");
320        # Start the table.
321        print SUMFILE "<table class=\"$tableStyle\">\n";
322        # Create the header row.
323        print SUMFILE Tr( { class => 'odd' }, th(["Group name",
324                                                 "Genomes",
325                                                 "Protein Encoding Genes (PEGs)",
326                                                 "Named genes in subsystems",            # s0
327                                                 "Named genes not in subsystems",        # n0
328                                                 "Hypothetical genes in subsystems",     # s1
329                                                 "Hypothetical genes not in subsystems", # n1
330                                                 "RNAs",
331                                                   ])) . "\n";
332        # Set up a flag for the odd-even styling.
333        my $rowFlag = 0;
334        # Put in the data rows.
335        for my $groupName (sort keys %summaries) {
336            my $group = $summaries{$groupName};
337            # Compute the link for the current group.
338            my $groupLink = a({ href => $groupData{$groupName}->[0] }, $groupName);
339            # Create the table row.
340            my $rowHtml = join("",
341                               td($groupLink),
342                               td({ class => $numStyle }, Tracer::CommaFormat($group->{genomes})),
343                               td({ class => $numStyle }, Tracer::CommaFormat($group->{pegs})),
344                               td({ class => $counterStyle }, [ map { Tracer::CommaFormat($group->{$_}) } @columnTypes ]),
345                               td({ class => $numStyle }, Tracer::CommaFormat($group->{RNAs})),
346                              );
347            print SUMFILE Tr( { class => $rowStyle[$rowFlag] }, $rowHtml ) . "\n";
348            # Flip the row style.
349            $rowFlag = 1 - $rowFlag;
350        }
351        # Terminate the table and close the file.
352        print SUMFILE "</table>\n";
353        close SUMFILE;
354      # We're all done.      # We're all done.
355      Trace("Processing complete.") if T(2);      Trace("Processing complete.") if T(2);
356  }  }
357    
 =head3 Fix  
   
 C<< my %fixedHash = Fix(%groupHash); >>  
   
 Prepare a genome group hash for processing. Groups with the same primary name will be combined.  
 The primary name is the first capitalized word in the group name.  
   
 =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;  
         }  
         # Append this group's genomes into the result hash.  
         Tracer::AddToListMap(\%retVal, $realGroupID, @{$groupHash{$groupID}});  
     }  
     # Return the result hash.  
     return %retVal;  
 }  
   
   
358  1;  1;

Legend:
Removed from v.1.15  
changed lines
  Added in v.1.28

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3