#!/usr/bin/perl -w

use strict;

####################################################
## With this script we can obtain the sequence of ## 
##     a region of the genome from a GFF file.    ##
####################################################

my @introngff = "";                 # an array to read every line

my $start = "";                     # start of the 100 pb sequence

my $stop = "";                      # stop of first exon-flank

my $start2 = "";                    # start of second exon-flank

my $stop2 = "";                     # end of the 100 pb sequece

my $up = "";                        # first half of the sequence (50 pb)

my $down = "";                      # second half of the sequence (50 pb)

my $seq = "";                       # the sequence!!

my $chr = "";                       # name of the chromosome to look for its sequence in the db

my $entrada = $ARGV[0];

my $llegint = 0;                    # to read every line

my $path = "/disc8/genomes/H.sapiens/golden_path_200405hg17/chromFa";   # location of the genome db

my $resultat = "";                        # file where we will print the results

#######################################################################################################

print "Tipe the name of the OUTPUT (name.fa):";      # tip the name of the output file (name.fa) 
                                                     # that would contain the nucleotide sequence 

$resultat = <STDIN>;     

open (FILE,$entrada) || die "Can't find $entrada\n";

open (OUT,">$resultat") || die "Couldn't open exonpairs.fa\n";

while (<FILE>) {

  chomp $_;

  @introngff = split(/\t/,$_);                     # stores the line into an array using tabs 

  $chr    = $introngff[0];
  $start  = $introngff[3] - 50;                    # define the start and the end of the sequence and
  $stop   = $introngff[3] - 1;                     # store them into variables
  $start2 = $introngff[4] + 1;
  $stop2  = $introngff[4] + 50;

  $up = `chr_subseq $path/$chr.fa $start $stop`;       # serch the first half of the sequence in the 
  chomp $up;                                           # genome db using the postions (start and stop)
  
  $down = `chr_subseq $path/$chr.fa $start2 $stop2`;   # the same, but with the second half

  $seq = $up.$down;                                    # concatenate $up and $down to obtain the hole seq


  print OUT ">".$introngff[8]."|".$introngff[0]."|".$introngff[3]."|".$introngff[4]."|".$introngff[6]."\n".$seq."\n";
                                        # print out the results into the output file with fasta format
                                        # giving a name to every seq
  $llegint = $llegint + 1;              # go to the next line

}


close(FILE);

close(OUT);
