#!/usr/local/bin/perl5.6.0 -w

=head1 NAME

hth2embl - Splice-friendly helix-turn-helix motif search using EMBOSS
helixturnhelix

=head1 SYNOPSIS

 hth2embl filename
 hth2embl -c 4 filename
 hth2embl -c 0 -t 11 filename
 hth2embl -d 3.5 -t 1 filename

=head1 DESCRIPTION

This program searches CDS features for helix-turn-helix motifs using
the EMBOSS program helixturnhelix. The CDS features should be supplied
in the form of an EMBL file. The features will be translated using the
Bacterial translation table by default. The -t switch may be used to
specify alternative tables using the standard numbering
system. Currently features with /partial or /pseudo qualifiers are
skipped because of the phase problem (at some point I will add support
for this by searching a 3-frame translation in such cases).

The -d switch sets the standard deviation cutoff for the
helixturnhelix search (see the EMBOSS documentation). The -c switch
will allow the colour code of the features produced to be changed
(default colour is 8). The output is an EMBL feature table sent to
STDOUT.

=head1 METHODS

None

=head1 AUTHOR

Keith James (kdj@sanger.ac.uk)

=head1 COPYRIGHT

Copyright (C) 2001 Genome Research Limited. All Rights Reserved.

=head1 DISCLAIMER

This software is provided "as is" without warranty of any kind. It may
be used, redistributed and/or modified under the same conditions as
Perl itself.

=cut

use strict;
use Getopt::Std;
use Bio::PSU::SeqFactory;
use Bio::PSU::FeatureFactory;

my $hth      = "helixturnhelix";
my $hth_opts = "-filter -auto -stdout";

my $usage = "Usage: hth2embl [-cdt] <file name>

\t-c Colour code of features  (default  8)
\t-d Standard deviation cutoff (default 2.5 SD)
\t-t Translation table number (default 11)

\tTables available are:

\t 1) Standard
\t 2) Vertebrate Mitochondrial
\t 3) Yeast Mitochondrial
\t11) Bacterial
";

my (%opts, $cutoff, $colour, $table);
getopts('c:d:ht:', \%opts);

defined $opts{h} and die $usage;
defined $opts{c} and $colour  = $opts{c};
defined $opts{d} and $colour  = $opts{d};
defined $opts{t} and $table   = $opts{t};

# Default colour is green
$colour = 8 unless defined $colour;

# Default translation table is bacterial (11)
$table ||= '11';

# Default SD cutoff is 2.5
$cutoff = 2.5 unless defined $cutoff;

my $file = shift;

unless (defined $file) { die $usage }

# Create EMBL stream and gather one sequence
my $es = Bio::PSU::SeqFactory->make(-file => "<$file");
my $e  = $es->next_seq;

# Handler for printing features as they appear, rather than printing a
# big feature table at the end
my $fth = Bio::PSU::IO::FTHandler->new;
my $ff = Bio::PSU::FeatureFactory->new;

F: foreach my $f ($e->features)
{
    # Discard non-CDS features
    next F unless $f->key eq 'CDS';

    # Discard any CDS with /pseudo or /partial qualifier
    next F if $f->qexists('pseudo|partial');

    # Translate with specified translation table
    my $aa = $f->translate($table);

    # Set a temporary file
    my $time    = time;
    my $hth_tmp = "hth.$time";

    # Apply SD cutoff
    $hth_opts .= " -minsd $cutoff";

    open(HTH, "| $hth $hth_opts > $hth_tmp") or
	die "Unable to open pipe to $hth: $!\n";

    # Create Fasta stream writer for translations
    my $fs = Bio::PSU::SeqFactory->make(-fh     => \*HTH,
                                        -format => 'fasta');
    $fs->write_seq($aa);
    close(HTH);

    my ($score, $sd, $seq, $start, $end);

    open(DAT, "<$hth_tmp") or die "Unable to open $hth_tmp: $!\n";
    while (<DAT>)
    {
        if (/No hits above/)
        {
            close(DAT);
            unlink($hth_tmp);
            next F;
        }

        if (/^Score\s+(\d+)\s+(\(\S+\s+SD\))/)
        {
            ($score, $sd) = ($1, $2)
        }

        if (/Sequence:\s+(\S+)/)
        {
            $seq = $1
        }

        if (/^\s+(\d+)\s+(\d+)/)
        {
            ($start, $end) = ($1, $2)
        }
    }

    close(DAT);
    unlink($hth_tmp);

    unless (defined $score &&
            defined $sd    &&
            defined $seq   &&
            defined $start &&
            defined $end)
    {
        die "Unable to parse HTH result for CDS:\n" .
            $fth->make_feature_block($f);
    }

    my $length = $end - $start + 1;

    my $hthf = $ff->annotate_translation($f, $start, $length);
    $hthf->qadd("note", "Predicted helix-turn-helix motif with score $score $sd at aa $start-$end, sequence $seq");
    $hthf->qadd("colour", $colour);

    # Print features as they are created using the handler
    print $fth->make_feature_block($hthf);
}
