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

View of /FigKernelScripts/recall_pegs.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (download) (as text) (annotate)
Thu Dec 8 21:10:49 2011 UTC (7 years, 11 months ago) by overbeek
Branch: MAIN
CVS Tags: rast_rel_2014_0912, rast_rel_2014_0729, mgrast_release_3_1_2, mgrast_version_3_2, mgrast_dev_12152011, HEAD
Changes since 1.1: +1 -0 lines
add blocked robot

use strict;
use Data::Dumper;
use Carp;

=head1 recall_pegs

Recall peg assignments using Kmer annotations, if the peg is not in a subsystem and there has been
no manual annotation has been made in the last year.

------

Example:

    echo 83333.1 | recall_pegs > proposed.3.column.assignments

would produce a 6-column table.  The first column three columns make
the output usable as input to conditionally_assign (i.e., PEG IDs, the
second the current annotation, and the third the proposed annotation.
The last three contain values reflecting the strength of the kmer hits: 

     # of nonoverlapping hits (5 or more)
     # of hits
     difference in # of hits between the best score and the next (minimum of 10)


------

=head2 Output Format

The output is a 3-column table that can be used as input to the command

    conditionally_assign

The columns are [PEG,Current-annotation,New-annotation]

=cut

my %robot;
$robot{'rapid_propogation'} = 1;
$robot{'automated_assignments'} = 1;
$robot{'pseed_update'} = 1;
$robot{'fig'} = 1;
$robot{'gdpusch'} = 1;
$robot{'Kmers'};

use FIG;
my $fig = new FIG;
use ANNOserver;
my $annoO = ANNOserver->new;

while (defined($_ = <STDIN>))
{
    if ($_ =~ /^(\d+\.\d+)/)
    {
	my $g = $1;
	my %do_not_change;
	$/ = "\n//\n";
	my @ann = map { chomp; my($peg,$ts,$user,@rest) = split(/\n/,$_); 
                        &ok($ts) ? [$peg, $ts,$user,join(" ",@rest)] : () }
	          `cat $FIG_Config::organisms/$g/annotations`;
	$/ = "\n";
	foreach my $tuple (@ann)
	{
	    my($peg,$ts,$user,$comment) = @$tuple;
	    if (($comment =~ /(set.*function)|(assign)/i) &&
		(! $robot{$user}))
	    {
		$do_not_change{$peg} = 1;
	    }
	}

	my @seqs;
	foreach my $peg ($fig->all_features($g,'peg'))
	{
	    if (! $do_not_change{$peg})
	    {
		my @tmp = $fig->peg_to_subsystems($peg);
		if (@tmp > 0)
		{
		    $do_not_change{$peg} = 1;
		}
		else
		{
		    push(@seqs,[$peg,'',$fig->get_translation($peg)]);
		}
	    }
	}
	next if (@seqs < 1);

	my $handle = $annoO->assign_function_to_prot( -input => \@seqs,
						      -assignToAll => 0,
						      -seqHitThreshold => 5,
						      -scoreThreshold => 10,
						      -kmer => 8 );

	while (my $hit = $handle->get_next)
	{
	    my($peg,$func,$otu,$sc,$non_ov_hits,$ov_hits,undef) = @$hit;
	    my $actual = $fig->function_of($peg);
	    if ($func && (&strip($actual) ne &strip($func)))
	    {
		if (($non_ov_hits >= 5) && 
		    (! (($actual =~ /chloroplast/i) && ($func !~ /chloroplast/i))) &&
		    (! (($actual =~ /mitochondri/i) && ($func !~ /mitochondri/i))))
		{
		    print join("\t",($peg,$actual,$func,$non_ov_hits,$ov_hits,$sc)),"\n";
		}
	    }
	}
    }
}

sub ok {
    my($ts) = @_;

    my $year_of_sec = 365 * 24 * 60 * 60;
    my $curr = time;
    return ($curr - $ts) < $year_of_sec;
}

sub strip {
    my($x) = @_;

    $x =~ s/\s*\#.*$//;
    return $x;
}

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3