#!/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 = 'TAG|TAA|TGA';

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 - $i),
							-alphabet => "dna");
	my @nt = split(//, $seq->seq);
	
	my $orfStart = -1;
	my $orfEnd = -1;
	my $strand = 1;
	
	for(my $b=0; $b<@nt; $b+=3){
		my $codon = $nt[$b].$nt[$b+1].$nt[$b+2];
		
		if($orfStart == -1 && $codon =~ /^$startCodon$/){
			$orfStart = $b+$i+1;
# 			print $codon."-".$startCodon."<BR>";
		}elsif($orfStart >=0 && $codon =~ /^$endCodon$/){
			$orfEnd = $b+3+$i;
			
			if ($orfEnd - $orfStart + 1 > $minSize){
# 				print $codon."-".$endCodon."*<BR>";
				push(@orfs, marcaOrf($ctgSeq, $orfStart, $orfEnd, $strand));
			}
			$orfStart = -1;
			$orfEnd = -1;
		}
	}
# 	exit;
}

# $endCodon = reverse $endCodon;
# $endCodon =~ tr/ATCG/TAGC/;
# 
# $startCodon = reverse $startCodon;
# $startCodon =~ tr/ATCG/TAGC/;

for(my $i=0; $i<3; $i++){
	my $seq = Bio::Seq->new(-seq => $ctgSeq->revcom->subseq($i+1, $ctgSeq->length),
				-alphabet => "dna");
	
# 	my $len = $revcom->length/3 + 1;
	
	my @nt = split(//, $seq->seq);
# 	for($i=0;$i<@nt;$i++){
# 		print $nt[$i]."\n";
# 	}
# 	exit;
	my $orfStart = -1;
	my $orfEnd = -1;
	my $strand = -1;
	for(my $b=0; $b<@nt; $b+=3){
		my $codon = $nt[$b].$nt[$b+1].$nt[$b+2];
		
		if($orfStart == -1 && $codon =~ /^$startCodon$/){
			print "i=$i b=$b $codon -".$startCodon."<BR>";
			$orfStart = $ctgSeq->length - $b - $i;
# $orfStart = $b+$i+1;
		}elsif($orfStart >=0 && $codon =~ /^$endCodon$/){
			$orfEnd = $ctgSeq->length - ($b+2) - $i;
# $orfEnd = $b+3+$i;
			if ($orfStart - $orfEnd + 1 > $minSize){
				print "<BR>".$codon." - ".$endCodon." ".$orfStart." ".$orfEnd."<BR>";
				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");
# 		$orfSeq = $orfSeq->revcom;
		$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);
	
	if($strand < 0){
# 	print $orfSeq->display_id."\n";
	print $orfSeq->seq."<BR>";
	}
	$orfHash{'seq'} = $orfSeq;

	
	return \%orfHash;
# 	exit;
}