#!/usr/bin/env perl

use strict;

my $usage = "$0 psl_file.psl [gff_output_file]

Psl files are generated by BLAT and have the following format (tab-separated):
n_match n_mismatch rep_match Ncount QgapCount QgapBases TgapCount TgapBases strand Qname Qsize Qstart Qend Tname Tsize Tstart Tend blockCount blockSizes qStarts tStarts
Last 3 fields are comma-separated lists. Nucleotides counted starting from 0.

Orf file is a tab-delimited file with the following columns
SEQ_name CDSstart CDSend CDSbreak + optional columns
CDSbreak is > 0 if CDS has a frame shift inside (e.g., due to errors in the genome or alignment errors)

To use psl2gif.pl with the NIA Gene Index, download files T-psl.txt and Torf-param.txt.
Then use command: psl2gif.pl blat_out.psl
";

my $arg=0;
my $input_file = $ARGV[$arg++] or die $usage;
my $output_file = 0;
$output_file = $ARGV[$arg++];

print "Executing: $0\n";
# my %hashChr;
# my $chr;
if (! $output_file){
    $output_file = $input_file;
    $output_file =~ s/\.psl$/\.gff/i;
}
#my $Torf_file = $input_file;
#$Torf_file =~ s/\.psl$/\.orf/i;

# open (INFO, $input_file) or die "file $input_file not found";
# while(my $line = <INFO>){
# 	my ($n_match,$n_mismatch,$rep_match,$Ncount,$QgapCount,$QgapBases,$TgapCount,$TgapBases,$strand,$Qname,$Qsize,$Qstart,$Qend,$Tname,$Tsize,$Tstart,$Tend,$blockCount,$blockSizes,$qStarts,$tStarts)=split(/\t/, $line);
# 	if(!$Tname){ next; }
# 	if(!$hashChr{$Tname}){
# 		$hashChr{$Tname}=1;
# 	}
# }
# close INFO;
# Get ORF information
my %hashORF;
open (INFO, $input_file) or die $!;
while(my $line = <INFO>){
    # my($Tname,$leftORF,$rightORF,$break,$ORFlen,$firstAA,$Kozak) = split (/\t/, $line);
    my @inf = split (/\t/, $line);
    $hashORF{$inf[9]} = [$inf[15],$inf[16],0];
}
close INFO;

open (OUTPUT, ">$output_file") or die "cannot create output file";
print "Writing output to $output_file\n";

#foreach my $chr (sort keys %hashChr){
# 	my $gff_file = "$chr.gff";
# 	if($output_file){
# 		$gff_file = $output_file . "_" . $gff_file;
# 	}
# 	open (OUTPUT, ">$gff_file") or die "cannot create output file";
# 	print "Writing output to $gff_file\n";
	open (INFO, $input_file) or die "file $input_file not found";
	while(my $line = <INFO>){
		my ($n_match,$n_mismatch,$rep_match,$Ncount,$QgapCount,$QgapBases,$TgapCount,$TgapBases,$strand,$Qname,$Qsize,$Qstart,$Qend,$Tname,$Tsize,$Tstart,$Tend,$blockCount,$blockSizes,$qStarts,$tStarts)=split(/\t/, $line);
		#if($chr ne $Tname){ next; }
		my ($leftORF,$rightORF,$break) = (-1,-1,-1);
		if($hashORF{$Qname}){
			($leftORF,$rightORF,$break) = @{$hashORF{$Qname}};
		}
		my @blockSize = split(/,/,$blockSizes);
		my @tStart = split(/,/, $tStarts);
		my @qStart = split(/,/, $qStarts);
		my $CDSstart = -1;
		my $CDSend = -1;
		my $frame_old;
		for(my $i=0; $i<@tStart; ++$i){
			my $start = $tStart[$i]+1;
			my $end = $tStart[$i]+$blockSize[$i];
			my $frameStart = $qStart[$i];
			my $qgap = 0;
			if($i>0){
				$qgap = $qStart[$i]-($qStart[$i-1]+$blockSize[$i-1]);
			}
			if($strand eq "-"){
				$frameStart = $Qsize-($qStart[$i]+$blockSize[$i]);
			}
			my $frameEnd = $frameStart+$blockSize[$i];
			my $frame = ".";
			if($frameStart >= $leftORF && $frameStart < $rightORF){
				$frame = ($leftORF-$frameStart)%3;
				if($break && $frameStart>=$break){
					$frame = ($break-$frameStart)%3;
				}
			}
			if($strand eq "+"){
				if($frameStart-$qgap <= $leftORF && $frameEnd > $leftORF){
					if($leftORF>=$frameStart){
						$CDSstart = $tStart[$i]+($leftORF-$frameStart)+1;
					}else{
						$CDSstart = $tStart[$i]+$frame+1;
					}
				}
				if($frameStart-$qgap < $rightORF && $frameEnd >= $rightORF){
					if($rightORF>$frameStart){
						$CDSend = $tStart[$i]+($rightORF-$frameStart);
					}else{
						$CDSend = $tStart[$i-1]+int(($blockSize[$i-1]-$frame_old)/3)*3+$frame_old;
					}
				}
			}else{
				if($frameStart <= $leftORF && $frameEnd+$qgap > $leftORF){
					if($frameEnd > $leftORF){
						$CDSend = $tStart[$i]+$blockSize[$i]-($leftORF-$frameStart);
					}else{
						$CDSend = $tStart[$i-1]+$blockSize[$i-1]-$frame_old;
					}
				}
				if($frameStart < $rightORF && $frameEnd+$qgap >= $rightORF){
					if($frameEnd >= $rightORF){
						$CDSstart = $tStart[$i]+$blockSize[$i]-($rightORF-$frameStart)+1;
					}else{
						$CDSstart = $tStart[$i]+$blockSize[$i]-(int(($blockSize[$i]-$frame)/3)*3+$frame)+1;
					}
				}
			}
			print OUTPUT "$Tname\tpsl2gff\texon\t$start\t$end\t.\t$strand\t$frame\t$Qname\n" if $end != 0;
			$frame_old = $frame;
		}
# 		if($CDSstart >= 0 && $CDSend >= 0){
# 			print OUTPUT "$Qname\tNIA_Mouse_4.0\tCDS\t$CDSstart\t$CDSend\t.\t$strand\t0\n";
# 		}
	}
	close INFO;
	# close OUTPUT;
#}
close OUTPUT;
exit(0);



