# -*- 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);