[Bio] / FigKernelScripts / call_genome_using_glimmer.pl Repository:
ViewVC logotype

View of /FigKernelScripts/call_genome_using_glimmer.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.3 - (download) (as text) (annotate)
Thu Jul 9 15:08:32 2009 UTC (10 years, 7 months ago) by gdpusch
Branch: MAIN
CVS Tags: mgrast_dev_08112011, mgrast_dev_08022011, rast_rel_2014_0912, myrast_rel40, mgrast_dev_05262011, mgrast_dev_04082011, rast_rel_2010_0928, mgrast_version_3_2, mgrast_dev_12152011, mgrast_dev_06072011, rast_rel_2009_0925, rast_rel_2010_0526, rast_rel_2014_0729, mgrast_dev_02212011, rast_rel_2010_1206, mgrast_release_3_0, mgrast_dev_03252011, rast_rel_2010_0118, rast_rel_2011_0119, mgrast_release_3_0_4, mgrast_release_3_0_2, mgrast_release_3_0_3, mgrast_release_3_0_1, mgrast_dev_03312011, mgrast_release_3_1_2, mgrast_release_3_1_1, mgrast_release_3_1_0, mgrast_dev_04132011, mgrast_dev_04012011, rast_rel_2009_07_09, rast_rel_2010_0827, myrast_33, rast_rel_2011_0928, mgrast_dev_04052011, mgrast_dev_02222011, mgrast_dev_10262011
Changes since 1.2: +51 -12 lines
Rewrite to remove overlap removal.

# -*- perl -*-
########################################################################
# Copyright (c) 2003-2006 University of Chicago and Fellowship
# for Interpretations of Genomes. All Rights Reserved.
#
# This file is part of the SEED Toolkit.
#
# The SEED Toolkit is free software. You can redistribute
# it and/or modify it under the terms of the SEED Toolkit
# Public License.
#
# You should have received a copy of the SEED Toolkit Public License
# along with this program; if not write to the University of Chicago
# at info@ci.uchicago.edu or the Fellowship for Interpretation of
# Genomes at veronika@thefig.info or download a copy from
# http://www.theseed.org/LICENSE.TXT.
########################################################################

use strict;
use warnings;

use FIG;
my $fig = FIG->new();

use NewGenome;
use GenomeMeta;
use Getopt::Long;

$0 =~ m/([^\/]+)$/o;
my $self  = $1;
my $usage = qq($self [--glimmerV=[2,3] --code=num --meta=metafile] OrgDir);

if (@ARGV && ($ARGV[0] =~ m/-h(elp)?/)) {
    print STDERR "\n   usage: $usage\n\n";
    exit (0);
}


my $code            = 11;
my $meta_file       = qq();
my $glimmer_version = 3;

my $rc   = GetOptions(
		      "glimmerV=s" => \$glimmer_version,
		      "code=s"     => \$code,
		      "meta=s"     => \$meta_file,
		      );
if (!$rc) {
    die "\n   usage: $usage\n\n";
}

$ENV{RAST_GLIMMER_VERSION} = $glimmer_version;

my $OrgDir = shift @ARGV;
(-d $OrgDir) || die "\n   $usage\n\n";

my $to_call  = NewGenome->new($OrgDir);
my $taxon_ID = $to_call->get_taxid();


#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#...Open auxiliary file
#-----------------------------------------------------------------------
my $called_by_file = qq($OrgDir/called_by);
open(CALLED_BY, ">>$called_by_file")
    || die "Could not append-open $called_by_file";


#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#...Ask the $to_call object to find the RNAs
#-----------------------------------------------------------------------
if (not $to_call->get_fids_for_type(qq(rna))) {
    my @tbl = $to_call->find_rna();
    print STDERR (qq(\n>>> Found ), (scalar @tbl), qq( probable RNA genes\n)) if $ENV{VERBOSE};
    
    foreach my $entry (@tbl) {
	chomp $entry;
	
	my ($fid, $loc, $func) = split(/\t/, $entry);
	
	$loc = $to_call->from_seed_loc($loc);
	if (not defined($fid = $to_call->add_feature( -type => 'rna',
						      -loc  => $loc,
						      -func => $func
						      )
			)
	    ) {
	    die qq(Could not add RNA feature: $loc, \'$func\'\n);
	}
    }
    
    $to_call->export_features();
}



#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#...Call genome with GLIMMER (no overlap cleanup)
#-----------------------------------------------------------------------

my $cmd  = qq();
if    ($glimmer_version == 2) {
    $cmd = "$FIG_Config::bin/run_glimmer2  $taxon_ID  $OrgDir/contigs  -code=$code";
}
elsif ($glimmer_version == 3) {
    $cmd = "$FIG_Config::bin/run_glimmer3  -code=$code  $taxon_ID  $OrgDir/contigs";
}
else {
    confess "ERROR: GLIMMER-$glimmer_version is not supported";
}

my $err_tmp = "$FIG_Config::temp/glimmer.$$.stderr";
open(GLIM, "$cmd 2> $err_tmp |") or confess "Error opening glimmer pipe $cmd: $!\n";

my @tmp_tbl;
while (<GLIM>) {
    push(@tmp_tbl, $_);
}

if (!close(GLIM)) {
    open(ERR, "<$err_tmp");
    while (<ERR>) {
	print STDERR $_;
    }
    close(ERR);
    confess "Glimmer pipe close failed \$?=$? \$!=$!\n";
}
    
foreach my $entry (@tmp_tbl) {
    chomp $entry;
    
    my ($fid, $loc, $func) = split(/\t/, $entry);
    
    $loc = $to_call->from_seed_loc($loc);
    if (defined($fid = $to_call->add_feature( -type => 'peg',
					      -loc  => $loc,
					      )
		)
	) {
	print CALLED_BY qq($fid\t$self (GLIMMER-$glimmer_version)\n);
    }
    else {
	die qq(Could not add RNA feature: $loc, \'$func\'\n);
    }
}

$to_call->export_features();

close(CALLED_BY);


exit(0);

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3