#!/usr/bin/perl
# Name: gBank2cox.pl
# Authour: Jack Min
# Date:  Dec 14, 2006
#
# updated on:
#	Feb 6, 06. updated to retrieve the phyllum information
# 
#  this program is used to retrieve the cox1 gene sequences (DNAs, cds, and proteins) from mitochondria 
# genome sequences.
# 
# Input: multiple gene bank records (GenBank format) and fasta sequence;
# Output: genomic gene sequence in fasta format,
#		cDNA fasta (only coding sequences),
#		protein sequence
#  ? exon_coordinates file;
#		list of gi and Accession number. (this list will be used to compare future downloading to retrieve newly deposited genes).
#
# Note: the output file will use the name of the input file, with extion "genomic.fasta, cds.fasta, protein.fasta" (gene bank sequecne in fasta formate);
#	"cds" (only coding dna sequence).

if($#ARGV != 1)
{
    print "\n";
    print "	usage: gBank2cox.pl gBank_file fasta_file\n\n";
    print "	please note:\n\n";
    print "	4 files will be produced after running the program.\n";
    print "	the output files will have the same name as input file, and with different extensions: \n";
    print "		'.gene.fasta' (gene genomic sequences in fasta format);\n";
    print "		'.cds.fasta' (coding dna sequences in fasta format)\n";
    print "		'.protein.fasta' (cox1 protein in fasta format)\n";
    print "		'.acc_gi' (accession list for coding sequences)\n";
    exit;
    }

 print "it is running ...\n";

 # initialize all variables
  my $input = $ARGV[0];
  my $input2 = $ARGV[1];
  my $geneFile= $input.".gene.fasta"; #gene bank sequences in fasta format
  my $cdsFile = $input.".cds.fasta"; #coding dna sequence file
  my $proFile =$input.".protein.fasta"; #protein file
  my $accFile = $input.".acc_gi"; #accession and gi list file generated from the input file
 
################# to save sequences and definition ############
open (INPUT, $input2) || die ("Cannot open file $input2\n\n");
$value ='';
$key = '';
%giFun = ();
%giSeq = ();
$count = 0;
$genomeSeq = '';
while (<INPUT>)
{
        $line = $_; #$_, default input variable
        chomp ($line);
	  $line =~ s/[\r\f\n]//;
        
	  if ($line =~ /^>/)
     {
		$count++;
		if ($count >1)
		{
			$giSeq{$key} = $genomeSeq;
			#print $key."\t".$genomeSeq."\n";
			$genomeSeq = '';
		}						
		@tmpA = split(/\|/, $line);
            $key = $tmpA[1];
		$pos = index ($line, " ");
            $fun = substr ($line, $pos+1);
		$giFun{$key} = $fun;
	}
	elsif ($line =~ /[AGCTagct]/)
	{
		$line =~ s/\s+//g; #replace any space
		$genomeSeq .= $line;
	}
}
$giSeq{$key} = $genomeSeq;
#print $key."\t".$genomeSeq."\n";
$genomeSeq = '';

close (INPUT);

####################
  my @in =(); #store input
  my $item='';
  my @accession =();
  my $coxFlag=0;
  my $gFlag=0;
  my $cdsFlag=0;
  my $gFlag=0;
  my $geneCount=0;
  my $geneFlag=0;
  my $joinSeq='';
  my $joinSeq2='';
  my $originFlag=0;
  my $endFlag=0;
  my $seq='';
  my $sequence='';
  my $gDNA='';
  my @fasta=();
  my @exon=();
  my $cds = 0;
  @all =(); #to store all entries of input

# Check for files exsistence
# Open file
$int = 0;
open(INPUT, $input) || die("Cannot open file $input!");

$offSet = tell(INPUT);
while ($record = get_next_record(INPUT))
{# while
	@recordA = split(/\n/, $record);
   	$int++;
	print $int."\n";
	$featureFlag=0;
	$coxFlag = 0;
	$geneFlag = 0;
	foreach $item (@recordA) #foreach start
	{
	   chomp ($item);
		$item =~ s/[\r\f\n]//;
#Start parsing .. important details
 #Accession #.
   #   print $item."\n";
	if ($item =~ /DEFINITION/)
		{
			@tmp = split (/\s+/, $item);
			$def = $tmp[1];
		}
		elsif ($item=~ /^ACCESSION\s+/)
      {
			@tmp = split (/\s+/,$item);
         $accNum = $tmp[1];
	 		#push (@all, $head."\n");
		}
		
		elsif ($item=~ /^VERSION\s+/)
      {
			@tmp = split (/GI\:/,$item);
         $gi= $tmp[1];
			$fun = $giFun{$gi};
			$head = "gi|".$gi."|".$accNum;
	 	}
		elsif ($item=~ /^SOURCE\s+/)
      {
			@tmp = split (/\s+/,$item);
      	$source = $tmp[1];
		}
		elsif ($item=~ /ORGANISM/)
      {
			@tmp = split (/  /,$item);
         $org = pop (@tmp);
			#push (@funA, $gi."\t".$accNum."\t".$org."\t".$fun."\n");
		}
		elsif ($item=~ /Eukaryota; Fungi;/)
      {
			@tmp = split (/\;/,$item);
         $phyllum = $tmp[2];
			#push (@funA, $gi."\t".$accNum."\t".$org."\t".$fun."\n");
		}
		
		elsif ($item =~ /\/strain=\"/)
		{
			@tmp = split (/[""]/, $item);
			$strain = $tmp[1];
		}	
		if ($item =~ /gene/ && $item =~ /\.\./)
		{
			$genRange = ''; 
			@tmp1 = split (/\s+/, $item);
			$genRange = pop (@tmp1); #print "gene range: ".$genRange."\n"; 
			#print "LINE 161".$item."\n";
			if ($item =~ /join\(/)
			{
				$geneFlag = 1;
			}
			next;
		}
		if ($item =~ /(\d+\))/ && $geneFlag == 1)
		{
			$genRange .= $item; 		
		}	
		if ($item=~ /\/gene="cox1"/ || $item=~ /\/gene="COX1"/ || $item=~ /\/gene="COI"/ )
	   {
        		$coxFlag=1; $geneFlag = 0;
			#	print "gene range ln163: ".$genRange."\n";
		}
#Parsing out exons.
#Check for CDS first...
      if( $item=~/CDS /)
      {
			$transFlag = 0;
			if ($coxFlag == 1)
			{
				if($item=~ /join/)
            {
               	$joinFlag=1;
               	$joinSeq=$item;
				#	print "Ln171: ".$joinSeq."\n";
						next;
            }
               if ($item !~ /join/ && $item=~/\.\./)
               {
                   		#$singleFlag=1;
                		$joinSeq=$item;
               		next;
               }
         }  
		}#if CDS
 # if exon list continues onto second line
 #  use join
      elsif ($item=~ /\.\./)
		{
			if ($coxFlag==1)
        	{
        		if ($joinFlag == 1)
           	{
      			$joinSeq .= $item;
            }
			}
      }  
#get rid of unwanted information
      if ($item =~ /translation=/ && $coxFlag == 1)
      {
           #print $gi."\t".$joinSeq."\n";
				$joinFlag=0;
				$transFlag = 1;
				@tmp = split (/=\"/, $item);
				$proSeq = $tmp[1];
				push (@exon, "gi|".$gi."|".$org."\t".$joinSeq);
				#print $gi."\t".$joinSeq."\n";
				$genRange =~ s/\s+//g;
				push (@grA, "gi|".$gi."|".$org."\t".$genRange."\n"); #@grA: gene range array
				next;
      }
    
		if ($item =~ /\w\"/ || $item =~ /(\s+)\"/)
		{
			if ($transFlag ==1)
			{
				$coxFlag = 0;
				$transFlag = 0;
				$proSeq .=$item;
				$proSeq =~ s/\"//;
				$proSeq =~ s/\s+//g;
				push (@proA, ">gi|".$gi."|".$org."\n".$proSeq."\n");
         			$item='';
         			$joinSeq='';
         			last;
			}
		}
		if ($item=~ /\w/)
      {
         if ($transFlag ==1)
         {
              $proSeq .= $item; 
         }
      }
      if ($item=~/\/\//)
      {
          	$endFlag=1;
				#push (@funA, $gi."\t".$accNum."\t".$org."\t".$strain."\t".$phyllum."\t".$fun."\n");
				#$strain = '';
      }
    
   }#end of for each
#########
push (@funA, $gi."\t".$accNum."\t".$org."\t".$strain."\t".$fun."\n");
$strain = '';
$offSet = tell(INPUT);
}#end of while
close (INPUT);
###########

$funFile = $input.".fun"; 
open (fun, ">$funFile");
print fun  @funA;
close (fun);

open (OUT, ">$proFile");
print OUT  @proA;
close (OUT);

$grFile = $input.".gene.range";
open (OUT, ">$grFile");
print OUT  @grA;
close (OUT);

#get out exon values
#Initialise all variables for exon parsing.

    $item='';
    @exonValues=();
    @exonValues1=();
    @exArray=();
    @exArray1=();
    $item1 = ();
    $item2=();

#parsing each exon stored in exon
foreach $item (@exon)
{  #1
	chomp ($item);
	@tmp = split (/\t/,$item);
	$head=@tmp[0];
	$join=@tmp[1];
	$join=~ s/\s//g; #remove space
	if ($join=~ /join/)
	{ #2
		@exArray= split(/[()]/,$join);
		if ($join=~ /complement/)
		{ #3
			$join= $exArray[2];
			# print $join;
			@exArray1 =split (/,/,$join);
			$int = scalar @exArray1;
			for($t= $int-1; $t>=0; $t--)
			{
				$item1= $exArray1[$t] ;
				chomp($item1);
				@tmp1=split (/\.\./,$item1);
				$line = $head."\t".$tmp1[1]. "\t".$tmp1[0]."\n";
				push( @exonValues, $line);
				#print @exonValues;
			}
			push(@exonValues,"\n");
		}#3
		else
		{
			$join = $exArray[1];
			@exArray1 = split (/,/,$join);
			foreach $item1 (@exArray1)
			{#6
				chomp($item1);
				@tmp2=split(/\.\./,$item1);
				$line= $head."\t".$tmp2[0]."\t".$tmp2[1] ."\n";
				push (@exonValues, $line);
			}#6
			push (@exonValues,"\n");
		} #5
	}
#Split out complement numbers
	else
	{
		if($join=~/complement/)
		{
			my @comAry =split (/[()]/, $join);
			$tmp = $comAry[1];
			@cAry =split (/\.\./,$tmp);
			$line=$head."\t".$cAry[1]."\t".$cAry[0]."\n";
			push (@exonValues, $line."\n");
		}
		else
		{
      	my @temp2 = split(/\.\./, $join);
      	$line = $head."\t".substr($temp2[0],3)."\t".$temp2[1]."\n";    push (@exonValues, $line."\n");
      }
	}
}

# print out all exon values

$exonFile = "out.exon";
 open (EXON, ">$exonFile" );
 print EXON (@exonValues);
 close (EXON);

# parsing out the cds sequence using the exons and genomic seq

$count = 0; #to indicate the first seq.
$flag2 = 0; #to indicate the 2nd gene, as some genomic seq has two genes.
$preId = ''; #the previous seq id
$id ='';
@cdsArray = (); #store cds
$cds =''; #each individual cds

$exonStart = 0;
$exonEnd = 0;

foreach  $exon (@exonValues)
{
	$count++;
	chomp ($exon);
	
	if ($exon !~ /\d+/)
	{
		next;
	}
        @tmp = split(/\t/,$exon);
        $id = $tmp[0]; #seq id
	if ($id =~ /\|/)
	{
		@tmp2 = split(/\|/, $id);
		$oriId = $id; #original id
		$id = $tmp2[1];
		$flag2 = 1;
	}

        $exonStart = $tmp[1]; #exon start
	#print "exon start:".$exonStart."\n";
	$exonEnd = $tmp[2]; #exon end

	#print "exon: ".$exon."\t"."preID:".$preId."\t"."ID:".$id."\n";
	if ($exonStart =~ /[><]/)
	{
		$exonStart = substr($exonStart, 1); #skip the first letter > or <;
	}
	if ($exonEnd =~ /[><]/)
	{
		$exonEnd = substr($exonEnd, 1); #skip the first letter < or >;
	}

	$seq = $giSeq{$id};

	if ($exonStart < $exonEnd)
	{
		$exonSeq = substr($seq, $exonStart-1, $exonEnd - $exonStart +1);
	#	print "exon Seq:".$exonSeq."\n";
	}
	else
	{
		$exonSeq = substr($seq, $exonEnd-1, $exonStart - $exonEnd +1);
		#print "exon end".$exonEnd."\t"."exon start:".$exonStart."\t"."exon seq:".$exonSeq."\n";

		$revSeq = reverse($exonSeq);
		$revSeq =~ tr/ACGTacgt/TGCAtgca/;
		#print "reverse complement:".$revSeq."\n";
		$exonSeq =$revSeq;
	}

	if ($flag2 == 0)
	{
		if (($id ne $preId) && ($count >1))#a new seqID
		{
			push (@cdsArray, ">".$oriId."\n".$cds."\n");
			# print "preID:".$preID."cds: ".$cds."\n";

			$preId = $id;
			$cds = $exonSeq;
		}
		else
		{
			$preId = $oriId;
			$cds .= $exonSeq;
		}
	}

	else
	{
		if (($oriId ne $preId) && ($count >1))#a new seqID
		{
		#print $preId."\n";	
		@idA = split (/\|/, $preId);
			$id = "gi|".$idA[1]."|".$idA[0]; 
			push (@cdsArray, ">".$preId."\n".$cds."\n");
		#	print "cds: ".$cds."\n";
		#	print "preID:".$preID."cds: ".$cds."\n";

			$preId = $oriId;
			$cds = $exonSeq;
			$flag2 = 0;
		}
		else
		{
			$preId = $oriId;
			$cds .= $exonSeq;
			$flag2 = 0;
		}
	}
}

#push last
push (@cdsArray, ">".$preId."\n".$cds."\n");
#print "cds: ".$cds."\n";
open(CDNA,">$cdsFile");
print CDNA (@cdsArray);
close (CDNA);

# open(IN,">$accFile");
# print IN (@accession);
# close (IN);

# open(IN,">$allFile");
# print IN (@all);
# close (IN);
######################### to extract whole gene seq #######
foreach  $gene (@grA)
{
	$comFlag = 0;
	chomp ($gene);
	@tmp = split(/\t/,$gene);
   $id = $tmp[0]; #seq id
	if ($id =~ /\|/)
	{
		@tmp2 = split(/\|/, $id);
		$oriId = $id; #original id
		$gi = $tmp2[1];
	}
	$range = $tmp[1];
	$range =~ s/[<>]//g;
	######### handling join range ############
	if ($range =~ /join/)
	{
		$range =~ s/join//;
	}
	if ($range =~ /complement/)
	{
		$comFlag = 1;
	}
	$range =~ s/[()]//;
	@tmp3 = split (/\.\./, $range);   
	$geneStart = $tmp3[0]; #gene start
	#print "gene start:".$geneStart."\n";
	$geneEnd = pop(@tmp3); #gene end

	if ($geneStart =~ /[><]/)
	{
		$geneStart = substr($geneStart, 1); #skip the first letter > or <;
	}
	if ($geneEnd =~ /[><]/)
	{
		$geneEnd = substr($geneEnd, 1); #skip the first letter < or >;
	}

	$seq = $giSeq{$gi};

	if ($comFlag ==0)
	{
		$geneSeq = substr($seq, $geneStart-1, $geneEnd - $geneStart +1);
	}
	else
	{
		$geneSeq = substr($seq, $geneEnd-1, $geneStart - $geneEnd +1);
		$revSeq = reverse($geneSeq);
		$revSeq =~ tr/ACGTacgt/TGCAtgca/;
		#print "reverse complement:".$revSeq."\n";
		$geneSeq =$revSeq;
	}
	$genLength = length ($geneSeq);
	push (@geneArray, ">gi|".$gi."\t".$genLength."\n".$geneSeq."\n");
}
open(OUT,">$geneFile");
print OUT (@geneArray);
close (OUT);


##########################
print "The work is done!!!\n";
exit;

sub get_next_record
{
        my($fh) = @_;
        use strict;
        my($offSet);
        my($record);
        my($seperator) = $/;
        $/ = "//";
        $record =<$fh>;
        $/ =$seperator;
        return $record;
}


