[Bio] / SubsystemExtension / multiple_subsystem_extension.pl Repository:
ViewVC logotype

View of /SubsystemExtension/multiple_subsystem_extension.pl

Parent Directory Parent Directory | Revision Log Revision Log

Revision 1.3 - (download) (as text) (annotate)
Thu May 4 04:12:24 2006 UTC (12 years, 11 months ago) by overbeek
Branch: MAIN
Changes since 1.2: +17 -7 lines
RAE: Adding automatic subsystems extension and queue to Heikos tool


use strict;
use warnings;

use Getopt::Std;
use Term::ReadKey;
use IO::Handle;
use Data::Dumper;
use FIG;
use Subsystem;
use SubsystemExtension::AEOS;
use SubsystemExtension::ValidationInterface;
use SubsystemExtension::SubsystemCandidateFactory::MultiSEED;
use SubsystemExtension::ExtensionConfig qw (TEMPDIR);

use constant DEBUG => 0;
no warnings qw(redefine);

#  this is necessary if the script is started via rsh(1)
#  otherwise you won't see any output until the first <RETURN>

sub usage {
    print "\t -= multiple_subsystem_extension =-\n\n";
    print "This script is used to extend subsystems (either one or a list in a file) to one ore more genomes.\n";
    print "If neither the parameter -t nor -m is defined all genomes not yet present in the subsystem (excluding environmental samples) will be used\n\n";
    print "usage:\n multiple_subsystem_extension -s <subsystem name> -f <subsystem file> [-t <taxon_id> -m <taxon_id file> -a] -d <validation_directory> -e <e-value cutoff: e.g. 1.0e-10> -r <number of organisms to search in> -c <score cutoff default=150 >\n\n";

# global variables
our($opt_s, $opt_t, $opt_f, $opt_d, $opt_m, $opt_a, $opt_e, $opt_c, $opt_r);


if ((!$opt_s) && (!$opt_f)) {
    print "ERROR: Can't start script : no subsystem/filename given!\n";
    exit 1;

my $fig = new FIG;
my @subsystems;

push @subsystems, $opt_s if ($opt_s);

## collect the subsystem names from a flat file
# each line holds one subsystem name

if ($opt_f && -e $opt_f) {
    open SSF, "<$opt_f";
    while (<SSF>) {
	push @subsystems, $_;
    close SSF;

my $html = [];

my @missing_genomes;

# add single or multiple taxon ids (parsed from file) to 
# the list of missing genomes

if ($opt_t) {
    push @missing_genomes, $opt_t;
} elsif ($opt_m && -e $opt_m) {
    open TF, "<$opt_m";
    while (<TF>) {
	push @missing_genomes, $_ if $_ =~ /\d+\.\d+/;
    close TF;
} elsif ($opt_a) {    # extend to all complete genomes
    @missing_genomes = grep { $_ !~ /^99999/ } $fig->genomes("complete");

$opt_d = TEMPDIR unless ($opt_d);

unless (-d $opt_d) {mkdir $opt_d}

foreach (@subsystems) {

    my $subsystem = $fig->get_subsystem($_);

    if (!$subsystem) {
        print STDERR "could not initialize the subsystem $opt_s\n";

    # RAE: if the directory to write the results in exists, rename it
    if (-d $opt_d."/".$subsystem->get_name())
        my $count=1;
        $count++ while (-e $opt_d."/".$subsystem->get_name().".".$count);
        rename $opt_d."/".$subsystem->get_name(), $opt_d."/".$subsystem->get_name().".".$count;
    mkdir $opt_d."/".$subsystem->get_name();
    my $extension = SubsystemExtension::SubsystemCandidateFactory::MultiSEED->new($subsystem, "master", $opt_e, $opt_r);

    if (ref $extension && $extension->isa("SubsystemExtension::SubsystemCandidateFactory")) {

		my %in = map { $_ => 1 } $subsystem->get_genomes;

		unless (scalar @missing_genomes > 0) {
			print STDERR "getting missing genomes for subsystem ".$subsystem->get_name()."\n" if (DEBUG); 
			@missing_genomes = grep { ! $in{$_} } grep { $_ !~ /^99999/ } $fig->genomes("complete");

		my $subsystem_candidates = $extension->generate(\@missing_genomes);
		foreach my $genome (keys %$subsystem_candidates) {
			my $subsystem_candidate;
		# generate the aeos and extend the subsystem!!
			my $aeos = SubsystemExtension::AEOS->new($fig, $subsystem, $genome, $subsystem_candidates->{$genome}, {mute=>1, score_cutoff => $opt_c});
			$subsystem_candidate = $aeos->extend();
			my $subsystem_name = $subsystem->get_name();
			$subsystem_name =~ s/\(/\\\(/;
			$subsystem_name =~ s/\)/\\\)/;
			$subsystem_candidate->toFile($opt_d."/".$subsystem->get_name()."/$genome.scs") if (-d $opt_d."/".$subsystem->get_name());

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3