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

View of /FigKernelScripts/make_PHOBs_for_subsystems.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.5 - (download) (as text) (annotate)
Mon Dec 5 18:56:37 2005 UTC (13 years, 11 months ago) by olson
Branch: MAIN
CVS Tags: mgrast_dev_08112011, rast_rel_2009_05_18, mgrast_dev_08022011, rast_rel_2014_0912, rast_rel_2008_06_18, myrast_rel40, rast_rel_2008_06_16, mgrast_dev_05262011, rast_rel_2008_12_18, mgrast_dev_04082011, rast_rel_2008_07_21, rast_rel_2010_0928, rast_2008_0924, mgrast_version_3_2, mgrast_dev_12152011, rast_rel_2008_04_23, mgrast_dev_06072011, rast_rel_2008_09_30, rast_rel_2009_0925, rast_rel_2010_0526, rast_rel_2014_0729, mgrast_dev_02212011, rast_rel_2010_1206, caBIG-05Apr06-00, mgrast_release_3_0, mgrast_dev_03252011, rast_rel_2010_0118, mgrast_rel_2008_0924, mgrast_rel_2008_1110_v2, rast_rel_2009_02_05, rast_rel_2011_0119, mgrast_rel_2008_0625, 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, rast_rel_2008_10_09, mgrast_dev_04012011, rast_release_2008_09_29, mgrast_rel_2008_0806, mgrast_rel_2008_0923, mgrast_rel_2008_0919, rast_rel_2009_07_09, rast_rel_2010_0827, mgrast_rel_2008_1110, myrast_33, rast_rel_2011_0928, rast_rel_2008_09_29, mgrast_rel_2008_0917, rast_rel_2008_10_29, mgrast_dev_04052011, mgrast_dev_02222011, caBIG-13Feb06-00, rast_rel_2009_03_26, mgrast_dev_10262011, rast_rel_2008_11_24, rast_rel_2008_08_07, HEAD
Changes since 1.4: +17 -0 lines
Add license words.

# -*- 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 FIG;
my $fig = new FIG;

$usage = "usage: make_PHOBs_for_subsystems Dir [SimilarityCutoff FracCoverage] < FileOf_SS-Role-PEG-Func";

(
 ($dir = shift @ARGV) && (! -d $dir)
)
    || die $usage;

if (@ARGV == 2)
{
    $sim_cutoff = shift @ARGV;
    $frac_cov   = shift @ARGV;
}

mkdir($dir,0777) || die "could not make $dir";
while (defined($_ = <STDIN>))
{
    chop;
    ($ss,$role,$peg) = split(/\t/,$_);
    push(@{$subH{$ss}->{$role}},$peg);
}

$tmpF = "/tmp/pegs$$.fasta";

open(INDEXSS,">$dir/index") || die "could not open $dir/index";
$ssN = 0;
foreach $ss (sort keys(%subH))
{
    ++$ssN;
    mkdir("$dir/$ssN",0777) || die "could not make $dir/$ssN";
    print INDEXSS "$ssN\t$ss\n";

    open(INDEXR,">$dir/$ssN/index") || die "could not open $dir/$ssN/index";
    $roleN = 0;
    foreach $role (sort keys(%{$subH{$ss}}))
    {
	++$roleN;
	
	$pegs = $subH{$ss}->{$role};
	if (@$pegs < 8)
	{
	    print STDERR "Skipping $dir/$ssN/$roleN, role $role --- only ", (scalar @$pegs), " pegs\n";
	    next;
	}
	
	print INDEXR "$roleN\t$role\n";
	
	open(TMP,">$tmpF") || die "could not open $tmpF";
	foreach $peg (@$pegs)
	{
	    if ($pseq = $fig->get_translation($peg))
	    {
		print TMP ">$peg\n$pseq\n";
	    }
	    else
	    {
		print STDERR "failed to get translation for $peg\n";
	    }
	}
	close(TMP);
	
	mkdir("$dir/$ssN/$roleN",0777) || die "could not make $dir/$ssN/$roleN";
	&FIG::run("split_and_trim_sequences $dir/$ssN/$roleN/split_info $sim_cutoff $frac_cov < $tmpF");
	
	if (-s "$dir/$ssN/$roleN/split_info/set.sizes")
	{
	    open(SZ,"<$dir/$ssN/$roleN/split_info/set.sizes") || die " could not open $dir/$ssN/$roleN/split_info/set.sizes";
	    while (defined($_ = <SZ>))
	    {
		if (($_ =~ /^(\d+)\t(\d+)/) && ($2 > 3))
		{
		    $n = $1;
		    &FIG::run("make_phob_from_seqs $dir/$ssN/$roleN/$n < $dir/$ssN/$roleN/split_info/$n");
		}
	    }
	    close(SZ);
	}
	else
	{
	    system("rmdir $dir/$ssN/$roleN") && warn "Could not rmdir $dir/$ssN/$roleN";
	    print STDERR "Skipping role $role in $dir/$ssN/$roleN\n";
	}
    }
    close(INDEXR);
}
close(INDEXSS);
unlink($tmpF);

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3