#!/usr/bin/perl
# longorf.pl v0208020920
# (c) Dan Kortschak 2002

use vars qw($USAGE);

use strict;
use Getopt::Long;
use Bio::SeqIO;

$USAGE = "longorf [--help] [--notstrict] [--verbose] [--graph] [--width printwidth] [--format seqformat] --input seqfile\n";

my ($sequencefile,$sequenceformat,$notstrict,$graph,$verb,$printwidth,$help) =
   (undef,        'fasta',        undef,      undef,undef,50,         undef);

&GetOptions('input|i=s'              => \$sequencefile,
            'format|f=s'             => \$sequenceformat,
            'notstrict|n'            => \$notstrict,
            'width|w=s'              => \$printwidth,
            'graph|g'                => \$graph,
            'verbose|v'              => \$verb,
            'help|h'                 => \$help,
            );

if ($help) {
   exec('perldoc', $0);
   die;
}

if (!defined $sequencefile) {
    die($USAGE . "\nPlease specify an input filename.\n");
}

sub longestORF {
   my $best=0;
   my ($bests,$beste,$beststrand)=(-1,-1,0);
   my $bestorf="";

   my $relaxed=$_[1];
   my $dna=Bio::Seq->new(-seq => $_[0]);
   my %strand=('+'=>$dna->seq,
               '-'=>$dna->revcom->seq);

   foreach my $direction (keys %strand) {
      my @starts=();
      my @ends=();
      if ($relaxed) {
         for (my $frame=0;$frame<3;$frame++) {
            unless ($strand{$direction}=~m/^.{$frame}(taa|tga|tag)/i) {
               push @starts,$frame+1;
            }
         }
      }
      while ($strand{$direction}=~m/(atg)/gi) {
         push @starts,pos($strand{$direction})-2;
      }

      while ($strand{$direction}=~m/(taa|tga|tag)/gi) {
         push @ends,pos($strand{$direction})-2;
      }
      push @ends,($dna->length-2,$dna->length-1,$dna->length);

      for my $s (@starts) {
         for my $e (@ends) {
            if ($e%3==$s%3 and $e>$s) {
               if ($e-$s>$best) {
                  $best=$e-$s;
                  ($bests,$beste,$beststrand)=($s,$e,$direction);
                  $bestorf=Bio::Seq->new(-seq=>$strand{$direction})->subseq($s,$e);
               }
               last
            } else {
               next
            }
         }
      }
   }
   return ($best,$bests,$beste,$beststrand,$bestorf);
}


my $seqio = new Bio::SeqIO('-format' => $sequenceformat,
                           '-file'   => $sequencefile );

my ($length,$start,$end,$direction,$sequence);
my $count=0;
my @lengths;
my $totallength=0;

while (my $dna = $seqio->next_seq) {
   $count++;
   ($length,$start,$end,$direction,$sequence)=longestORF($dna->seq,$notstrict);
   if ($verb) {
      print $dna->display_id," ",$dna->desc,": "; 
      print "$length, $start, $end ($direction)\n$sequence\n\n",Bio::Seq->new(-seq=>$sequence)->translate->seq,"\n\n--\n\n";
   }
   $totallength+=$length;
   $lengths[$length/3]++;
}

print "Average ORF length: ", $totallength/$count,"\n\n";

print "Length distribution is:\n";

if ($graph) {
   my $length;
   my $maxlength=0;
   for ($length=0;$length<@lengths;$length++) {
      $lengths[$length]=0 unless $lengths[$length];
      $maxlength=$lengths[$length] if ($lengths[$length]>$maxlength);
   }
   for ($length=0;$length<@lengths;$length++) {
      print $length*3,"\t",$lengths[$length],"\t|";
      print "#"x(($lengths[$length])*$printwidth/$maxlength);
      print "\n";
   }
} else {
   for ($length=0;$length<@lengths;$length++) {
      print $length*3,"\t",($lengths[$length]or"0"),"\n";
   }
}

__END__

=head1 NAME

longorf.pl - perl script to find the longest ORF of a sequence

=head1 SYNOPSIS

% longorf.pl [-h] [-n] [-v] [-g] [-w printwidth] [-f seqformat] -i seqfile

=head1 DESCRIPTION

This script will examine a set of nucleotide sequences and determine
the longest ORF in each sequence. ORFs may start at the canonical ATG
or at the beginning of the sequence if the notstrict option is chosen.
The script will output a list of the longest ORF lengths, starts, ends
and strands with the ORF and amino acid sequence if the verbose option
is chosen. A histogram of the longest ORFs in the input set may be
printed by choosing the graph option.

=head1 FEEDBACK

This script is not supported by anyone, but requests can be made to the
author.

=head1 AUTHOR - Dan Kortschak <kortschak@rsbs.anu.spanner.edu.au>

=cut