#!/usr/bin/perl -w

use strict;
use DBI;
use DBD::mysql;
use Bio::Seq;
use Bio::SeqIO;
use Data::Dumper;

## MARCA TODAS AS ORFs DO CONTIG
## salva as ORF no arquivo genes.fas estilo Glimmer
## >orf00002  46 1404  len=1359

my $ctgFileName = shift;
my $startCodon = shift;
my $minSize = shift;

my $ctgSeq = Bio::SeqIO->new(-file => "<".$ctgFileName,
     	                     -format => "fasta");

$ctgSeq = $ctgSeq->next_seq();

my @seqForwardFrames;

my $endCodon = '\*';

my @orfs;

open (PREDICT, '>genes.predict.tmp');
print PREDICT ">contig\n";
close (PREDICT);  

for(my $i=0; $i<3; $i++){
	my $seq = Bio::Seq->new(-seq => $ctgSeq->subseq($i+1, $ctgSeq->length),
							-alphabet => "dna");
	
	my $translation = $seq->translate->seq;
# 	print $translation; exit;
	@{$seqForwardFrames[$i]} = split(//, $translation);
# 	print $seqForwardFrames[$i][0]."\n";
	
	my @aa = split(//, $translation);
	my $orfStart = -1;
	my $orfEnd = -1;
	my $strand = 1;
	for(my $b=0; $b<@aa; $b++){
		if($orfStart == -1 && $aa[$b] =~ /^$startCodon$/){
			$orfStart = $b*3 + $i + 1;
# 			print $aa[$b]."\n";
		}elsif($orfStart >=0 && $aa[$b] =~ /^$endCodon$/){
			$orfEnd = ($b+1)*3 + $i;
			if ($orfEnd - $orfStart + 1 > $minSize){
				push(@orfs, marcaOrf($ctgSeq, $orfStart, $orfEnd, $strand));
			}
			$orfStart = -1;
			$orfEnd = -1;
		}
	}
# 	exit;
}

for(my $i=0; $i<3; $i++){
	my $seq = Bio::Seq->new(-seq => $ctgSeq->revcom->subseq($i+1, $ctgSeq->length),
				-alphabet => "dna");
	
	my $translation = $seq->translate->seq;
# 	print $translation; exit;
	@{$seqForwardFrames[$i]} = split(//, $translation);
# 	print $seqForwardFrames[$i][0]."\n";
	
	my @aa = split(//, $translation);
	my $orfStart = -1;
	my $orfEnd = -1;
	my $strand = -1;
	for(my $b=0; $b<@aa; $b++){
		if($orfStart == -1 && $aa[$b] =~ /^$startCodon$/){
			$orfStart = $ctgSeq->length - $b*3 - $i;
		}elsif($orfStart >=0 && $aa[$b] =~ /^$endCodon$/){
			$orfEnd = $ctgSeq->length - ($b+1)*3 - $i + 1;
			if ($orfStart - $orfEnd + 1 > $minSize){
				push(@orfs, marcaOrf($ctgSeq, $orfStart, $orfEnd, $strand));
			}
			$orfStart = -1;
			$orfEnd = -1;
		}
	}
# 	exit;
}

system("rm genes.fas");

@orfs = sort {$a->{'start'} <=> $b->{'start'}} @orfs;

foreach my $orf (@orfs){
	my $seq_out = Bio::SeqIO->new(-format => 'fasta',
				   -file => ">>genes.fas");

	$seq_out->write_seq($orf->{'seq'});   
}

system("cat genes.predict.tmp | sort -n -k2 > genes.predict");

sub marcaOrf{
	my ($ctgSeq, $start, $end, $strand) = @_;
	
	my $orfSeq;
	my $orfId;
	my %orfHash;
	
	$orfHash{'start'} = $start;

	if ($strand > 0){
		$orfSeq = Bio::Seq->new(-seq => $ctgSeq->subseq($start,$end),
				   -alphabet => "dna");

		$orfId = "orf $start $end  len=".($end-$start+1);
		$orfSeq->display_id("$orfId");
	}else{
		$orfSeq = Bio::Seq->new(-seq => $ctgSeq->subseq($end,$start),
				   -alphabet => "dna");
		
		$orfId = "orf $start $end  len=".($start-$end+1);
		$orfSeq->display_id("$orfId");
	}
	
	open (PREDICT, '>>genes.predict.tmp');
	print PREDICT "orf00000    $start    $end\n";
	close (PREDICT);
	
# 	print $orfSeq->display_id."\n";
# 	print $orfSeq->seq."\n";
	$orfHash{'seq'} = $orfSeq;
	
	return \%orfHash;
# 	exit;
}