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

View of /FigKernelScripts/blocked_by.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (download) (as text) (annotate)
Mon Feb 24 21:38:50 2014 UTC (5 years, 8 months ago) by overbeek
Branch: MAIN
CVS Tags: rast_rel_2014_0729, rast_rel_2014_0912, HEAD
find out why a kmer was blocked

#
# 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 Data::Dumper;

use strict;
use SeedEnv;
use FIG_Config;

use FIG;
my $fig = new FIG;

# usage: blocked_by PEG > potenitally conflicting PEGs/annotations

my $peg1 = $ARGV[0];
$peg1 || die "give a PEG";

my $func1 = &SeedUtils::strip_func_comment(scalar $fig->function_of($peg1));
my $genome = &SeedUtils::genome_of($peg1);
my $seq1   = $fig->get_translation($peg1);
my %kmers_in_peg1;
my %kmers_in_peg1_hit;
for (my $i=0; ($i < (length($seq1)-8)); $i++)
{
    my $kmer = uc substr($seq1,$i,8);
    $kmers_in_peg1{$kmer} = 1;
}

my %hits;
# my @genomes = map { (($_ =~ /(\d+\.\d+)$/) && ($1 != $genome)) ? $1 : () } `echo $genome | svr_members_of_otu`;
my @genomes = grep { $fig->is_prokaryotic($_) } $fig->genomes('complete');
foreach my $g (@genomes) 
{
    if (-d "$FIG_Config::organisms/$g")
    {
	print STDERR "processing $g\n";
	my @pegs = $fig->all_features($g,'peg');
	foreach my $peg2 (@pegs)
	{
	    my $func2 = &SeedUtils::strip_func_comment(scalar $fig->function_of($peg2));
	    if ($func2 ne $func1)
	    {
		my $seq2   = $fig->get_translation($peg2);
		my $ln = length($seq2);
		for (my $i=0; ($i < ($ln-8)); $i++)
		{
		    my $kmer = uc substr($seq2,$i,8);
		    if ($kmers_in_peg1{$kmer})
		    {
			print STDERR "hit: $peg2\n";
			$kmers_in_peg1_hit{$kmer} = 1;
			$hits{$peg2}++;
		    }
		}
	    }
	}
    }
}

foreach my $peg ( sort { $hits{$b} <=> $hits{$a} } keys(%hits))
{
    print join("\t",($hits{$peg},$peg,scalar $fig->function_of($peg))),"\n";
}
print "====\n";
foreach my $kmer (sort keys(%kmers_in_peg1))
{
    if (! $kmers_in_peg1_hit{$kmer})
    {
	print $kmer,"\n";
    }
}

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3