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

=head1 NAME

tmhmm2embl - Generate EMBL annotation of transmembrane helices
predicted by TMHMM2.0

=head1 SYNOPSIS

Examples:

 tmhmm2embl -t 11 embl_file > tmhmm.tab
 tmhmm2embl -lv tmhmm_file > tmhmm.tab
 tmhmm2embl -c 4 -t 1 embl_file > tmhmm.tab

The EMBL feature table is written to STDOUT.

=head1 DESCRIPTION

This program requires an EMBL file containing annotated CDS
features. The features will be translated using the specified
translation table (default is table 11; bacterial), except for those
with /pseudo, /partial or /transl_except qualifiers. TMHMM2.0 will be
run on these to produce at output file with a '.tmhmm' extension.

The output file is parsed by the script to produce an EMBL table on
STDOUT. It may also be saved and re-loaded to produce other tables by
running the tmhmm2embl with the -l switch.

At the moment gene names are used to identify each CDS. If no gene
name is found the token 'Unknown_CDS' will appear instead.

The -c switch sets the colour of the features (using the Artemis
colour codes) and the -v switch will provide verbose annotation of all
the helices.

The -s switch will create a separate feature for each predicted
transmembrane helix.

This script requires the Bio::PSU Perl modules see:

http://www.sanger.ac.uk/Users/kdj/software.html

Don't forget to cite Anders Krogh if you use this:

TMHMM is described in

Erik L.L. Sonnhammer, Gunnar von Heijne, and Anders Krogh:
A hidden Markov model for predicting transmembrane helices in protein
sequences.
In J. Glasgow et al., eds.: Proc. Sixth Int. Conf. on Intelligent
Systems for Molecular Biology, pages 175-182. AAAI Press, 1998.

=head1 METHODS

None

=head1 AUTHOR

Keith James (kdj@sanger.ac.uk)

=head1 COPYRIGHT

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

decodeanhmm, version 1.1e Copyright (C) 1998-1999 by Anders Krogh. 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;

# Configuration of TMHMM variables
my $bindir = "/nfs/disk222/yeastpub/bio-soft/TMHMM2.0/bin";
my $libdir = "/nfs/disk222/yeastpub/bio-soft/TMHMM2.0/lib";

my $tmhmm_bin = "$bindir/decodeanhmm";
my $tmhmm_mod = "$libdir/TMHMM2.0.model";
my $tmhmm_opt = "$libdir/TMHMM2.0.options";

# TMHMM execution args
my $tmhmm_args = "$tmhmm_mod -f $tmhmm_opt";

my (%opts, $colour, $load, $split, $table, $verbose);
my $usage = "Usage: tmhmm2embl [-l | -t <table>] [-c <colour> -v] file

\t-t specify translation table (standard numbering, incompatible with -l)
\t-l use pre-run TMHMM file instead of EMBL file (no translation required)
\t-s split result into a feature for each helix within a CDS
\t-c specify which Artemis colour code to use (default 0)
\t-v verbose annotation (reports coordinates of all helices)
";

# Process options
getopts('c:hlst:v', \%opts);

defined $opts{h} and die $usage;
defined $opts{c} and $colour  = $opts{c};
defined $opts{l} and $load    = $opts{l};
defined $opts{s} and $split   = $opts{s};
defined $opts{t} and $table   = $opts{t};
defined $opts{v} and $verbose = $opts{v};

# -l and -t flags are incompatible
if ($load and $table){ die $usage }

# Default colour is white
$colour = 0 unless defined $colour;

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

my $file = shift;

unless (defined $file) { die $usage }

if (! $load)
{
    my $embl = $file;

    # Read EMBL stream and get entry
    my $es = Bio::PSU::SeqFactory->make(-file => "<$embl");
    my $e  = $es->next_seq;

    # Create Fasta output stream
    my $fo = Bio::PSU::SeqFactory->make(-fh     => \*TMHMM,
										-format => 'fasta');

    # Open a pipe to TMHMM
    open(TMHMM, "| $tmhmm_bin $tmhmm_args > $embl.tmhmm 2>/dev/null") or die "Unable to open pipe to $tmhmm_bin: $!\n";

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

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

		# Find and fix the id
		my $id = $f->gene;

		$id ||= "Unknown_CDS";
		$id   =~ s/[,;:]$//;

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

		# Translation exceptions need to be skipped for now
		if (grep /\*/, $aa->str)
		{
			if ($f->qexists('transl_except'))
			{
				next
			}
			else
			{
				warn "Stop codon present in ID $id. Skipping this CDS\n";
				next;
			}
		}

		# Set the id and desc for the Fasta header
		$aa->id($id);
		$aa->desc(join(':', $f->start, $f->end, $f->strand));

		# Write translation to TMHMM pipe as Fasta
		$fo->write_seq($aa);
	}
    close(TMHMM);
}

if (! $load)
{
    # Create an EMBL feature table from the TMHMM output
    open(IN, "<$file.tmhmm") or die "Unable to open $file.tmhmm: $!\n";
}
else
{
    # Open the pre-run data file
    open(IN, "<$file") or die "Unable to open $file: $!\n";
}

my $fth = Bio::PSU::IO::FTHandler->new;

LINE: while (my $line = <IN>)
{
    if ($line =~ /^>(\S+) (?:\S+ )*(\S+)/)
    {
		my ($id, $desc) = ($1, $2);
		my ($start, $end, $strand) = split(':', $desc);

		my (@prediction, @ranges, @coords);

    HIT: while (my $hit = <IN>)
		{
			if ($hit =~ /^>/)
	    {
			die "Malformed TMHMM output file: no prediction for ID $id\n"
	    }

			if ($hit =~ /^\%pred NB\(0\):(.*)/)
			{
				@prediction = split(',', $1);
				last HIT;
			}
		}

    SPAN: foreach my $span (@prediction)
		{
			my ($type, $from, $to);

			if ($span =~ /([ioMO]) (\d+) (\d+)/)
			{
				($type, $from, $to) = ($1, $2, $3);

				push(@coords, "$from-$to") if $type eq 'M';
			}
			else
			{
				die "Unable to parse span data from [$span] for ID $id\n"
			}

			# Skip unless a predicted membrane span
			next SPAN unless $type eq 'M';

			unless ($start && $end && $strand)
			{
				die "Problem at $id: start, end or strand undefined\n"
			}

			# Calculate new start and end
			my ($nstart, $nend);
			if ($strand == 1)
			{
				$nstart = $start + (3 * ($from - 1));
				$nend   = $start + (3 * $to)   - 1;
			}
			else
			{
				$nstart = $end - (3 * $to)   + 1;
				$nend   = $end - (3 * $from);
			}

			my $range = Bio::PSU::Range->new(-start  => $nstart,
											 -end    => $nend,
											 -strand => $strand);
			push(@ranges, $range);
		}

		# Skip unless spans were found
		next LINE unless @ranges;

		if ($split)
		{
			foreach (0..$#ranges)
			{
				my $a = Bio::PSU::Feature->new(-key    => 'misc_feature',
											   -ranges => $ranges[$_]);

				my $note = "probable transmembrane helix predicted for $id by TMHMM2.0 at aa ". $coords[$_];
				$a->qadd('colour', $colour);
				$a->qadd('note', $note);

				# Print the feature, using order instead of join
				my $text = $fth->make_feature_block($a);
				$text =~ s/join/order/g;
				print $text;
			}
		}
		else
		{
			# Create new feature(s) with some annotation
			my $a = Bio::PSU::Feature->new(-key    => 'misc_feature',
										   -ranges => \@ranges);

			my $helixcount = scalar @ranges;
			my $helix      = $helixcount > 1 ? "helices" : "helix";

			my $note = "$helixcount probable transmembrane $helix predicted for $id by TMHMM2.0";

			if ($verbose)
			{
				$note .= " at aa ";

				my $text = join(', ', @coords);
				$text =~ s/, (\d+-\d+)$/ and $1/;

				$note .= $text;
			}

			$a->qadd('colour', $colour);
			$a->qadd('note', $note);

			# Print the feature, using order instead of join
			my $text = $fth->make_feature_block($a);
			$text =~ s/join/order/g;
			print $text;
		}
    }
}

close(IN);







