#!/usr/bin/perl -w

if ($#ARGV < 2) {die "ERROR need 3 parameters: \n".
		     "     genome input file (../genomes/NC_000913) in genbank format,\n".
		     "     coding indices (../genomes/NC_000913.coding in .ptt format)".
		     "     output-directory (../geblastet/blast_in/"; }

$MIN_LENGTH = 20;
$seq_file = $ARGV[0];
$seq_coding = $ARGV[1];
@tmp = split(/\//,$seq_file);
$seq_name = $tmp[$#tmp];
$org = substr($seq_name,0,(length $seq_name) - index($seq_name,".") + 1); # orgname ohne punkte
$seq_non_cd = $ARGV[2].$seq_name;
print "$seq_name $seq_file\n"; 
$file_name = $seq_name;
$genomelaenge = `grep 'LOCUS'  $seq_file \| awk '{print \$3;}'`;
print "output ".$ARGV[2].$org." \n"; 

# read genome file to @all :
open(FILE,"<".$seq_file) or die " file ".$seq_file." not found.\n";
print "- file ".$seq_file." open...";
@all = <FILE>;
close(FILE);
print ", read and closed.\n";

# find rna sequence in file:
$i=0;
do {
	 $line = $all[$i];
    $i = $i + 1;
} until (index($line,"ORIGIN") >= 0);


# save sequence to @seq
for (; $i < ($#all - 1); $i++) {
    $line = $all[$i];
    $line =~ s/[0-9]//g;
    $line =~ s/ //g;
    push(@seq,$line);
}


# find non-coding regions :
print "- cutting non-coding region...\n";

# read indices of protein-coding genes from *.coding
open(FILE,"<".$seq_coding) or die "couldn't find *.coding file \n";
@all = <FILE>;
close(FILE);

$i=0;

do {
    $line = $all[$i];
    $i = $i + 1;
} until (index($line,"Location") >= 0);


for (; $i <= $#all; $i++) {
    $line = $all[$i];
    $line =~ s/ //g;
    @begin = split(/\./,$line);
    @end=split(/[\+\-]/,$begin[2]);
    push(@indx,$begin[0]);
    push(@indx,$end[0]);
}

# sort indices:
# (sorting a hash)
print " ,sorting region indices, ... ";
for ($i=0; $i<$#indx; $i=$i+2) {
    $tmp{$indx[$i]} = $indx[$i+1];
}
@keys = sort { $a <=> $b } (keys %tmp);
$i = 0;
foreach $key (@keys) {
    $indx[$i] = $key;
    $i++;
    $indx[$i] = $tmp{$key};
    $i++;
}

# write $seq to list @allletters:
foreach(@seq) {
    $_ =~ s/\s//g;
    push(@allletters,split(//,$_));
}

# write non-coding regions to fasta-file:
open(FILE,">", $ARGV[2].$org);

$begin = 0;
$end = $indx[0]-2;
$id = 0;
$laenge = $end - $begin + 1;
if ($begin<=$end && $laenge > $MIN_LENGTH) {
    $line="";
    print FILE ">".$org.":0:".($begin+1).":+:".$genomelaenge;
	 $id++;
    for ($i=$begin;$i<=$end;$i++) {
	$line=uc($line.$allletters[$i]);
    }
    print FILE $line."\n";
}

$indxpos=-1;
do {
    $indxpos += 2;
    $begin = $indx[$indxpos];
    $end = $indx[$indxpos+1]-2;
	 $laenge = $end - $begin + 1;
    if ($begin<=$end && $laenge > $MIN_LENGTH) {

	print FILE ">".$org.":0:".($begin+1).":+:".$genomelaenge;
	$id++;
	$line="";
#    print $begin.",".$indxpos.",".$end.".";
	for ($i=$begin;$i<=$end;$i++) {
	    $line=uc($line.$allletters[$i]);
	}
	print FILE $line."\n";
    }
} while($indxpos < ($#indx-2));

$begin = $indx[$#indx];
$end = $#allletters;
$laenge = $end - $begin + 1;
if ($begin<=$end && $laenge > $MIN_LENGTH) {
    $line="";
    print FILE ">".$org.":0:".($begin+1).":+:".$genomelaenge;
	 $id++;
    for ($i=$begin; $i<=$end; $i++) {
	$line=uc($line.$allletters[$i]);
    }
    print FILE $line;
}
close(FILE);

print "done successfuly\n";
print "- multiple fasta file with non-coding regions is:\n  ".$seq_non_cd."\n";
