#!/usr/bin/perl -w

# M.chandrasekhar, Dept of. EE, IIT Bombay,
# chandra_m_sekhar@hotmail.com
# Created on 4th February 2006
#
# This searches ATG positions in each sequence having 100 upstream and 
# 100 downstream and dumps the 203 length ATG window sequences
# additionally the <sequence number,serial number of ATG(203 length window), its position>
# is dumped on STDOUT. The position 500 (array index starts with index 0) 
# in each full length sequence is positive TIS
# 
# Usage: perl getvalidsequences.pl <input sequence file> <output 203 window sequences file> >  mapfile
# e.g.: perl getvalidsequences.pl gametsequence.txt gametATGsequences.txt > mapgamet
use strict;
$|++;


my $numseq=0;
my @dnaseqs=();
my $seqlength=undef;
my @NCLASS=();
my $winlen=100;

sub searchstring
{
       my $tstring=shift;
       my $qstring=shift;
       my $idx=0;
       my $i=0;
       my $tidx;
       my @indexs;
       $indexs[$i]=-1;
       
       while( ($tidx = index($tstring,$qstring,$idx)) != -1)
       {
          $indexs[$i] = $tidx;
          $idx = $tidx + length($qstring);
          $i++;
       }

       return($i,\@indexs);
}
   

sub getvalidindices
{
       my $seq=shift;
       my $str2 = shift;

       my $seqlength1 = length($str2);
       my @valindx;
       $valindx[0]=-1;
       foreach (@$seq)
       {
          # for 100
          # if( (($_ + 102) <= ($seqlength - 1)) && (($_ - 100) >= 0) ) 
          # for 50
          #if( (($_ + 52) <= ($seqlength - 1)) && (($_ - 50) >= 0) )
          if( (($_ + $winlen + 2) <= ($seqlength1 - 1)) && (($_ - $winlen) >= 0) ) 
          {
             if($valindx[0] == -1)
             {
               $valindx[0] = $_;
             }else{
               push(@valindx,$_);
             }
          }
       }
       return(\@valindx);
}

sub loadFile
{
    my $filename=shift;
    my $count=1;
    my $prev_seqlength = undef;
    open(FHANDLE,"<$filename") || die "Error: opening file $filename $!";
    while(<FHANDLE>)
    {
       chomp;
       my $fstring=$_;
       $seqlength = length($fstring) if(not defined $seqlength);
       $prev_seqlength = length($fstring);
       
       #if($prev_seqlength != $seqlength){
          #die "Error Mismatch sequence length between seqid $count and ". ($count - 1);
          #}
       push(@dnaseqs,$fstring);
       $count++;
    }
    close(FHANDLE);
}

sub findindices
{
   my $qqstr=shift;
   my @seqindices=();
   foreach(@dnaseqs)
   {

       my ($numc,$refindexs) = &searchstring($_,$qqstr);
       push(@seqindices,$refindexs);
   }
   return(\@seqindices);
}


sub constructfeatures
{
   my $seq1 = shift;
   my @mapseq=();
   my $flankwin = 0;


   my $di=0;
   $flankwin = ($winlen * 2) + 3;
   foreach my $seq (@$seq1)
   {
      my $dstr = $dnaseqs[$di];
      foreach (@$seq)
      {
         if($_ != -1)
         {
            #my $winseq = substr($dstr,$_ - $winlen,203); # -100,ATG,+100
            my $winseq = substr($dstr,$_ - $winlen,$flankwin); # -100,ATG,+100
            push(@NCLASS,$winseq);
            my @xim=($di,$_);
            push(@mapseq,\@xim);
         }
      }
      $di++;
   }

   return(\@mapseq);
}
    
    die("error: supply the fasta file\n") if($#ARGV+1 != 2);
    &loadFile($ARGV[0]);
    my $nseq = scalar @dnaseqs;
    my $seqindices = &findindices("ATG");
    my @seq_indices = ();
   
    my $blu = 0;
    foreach my $seq (@$seqindices)
    {
       my $validindices = &getvalidindices($seq,$dnaseqs[$blu]);
       push(@seq_indices,$validindices);
       $blu++;
    }

    die "error mismatch in number of sequences" if ($nseq != scalar @seq_indices);

    my $map_seq = &constructfeatures(\@seq_indices);
    my @mapSeq = @$map_seq; # holds info in seqs in negative class to orignial seqs

    my $outfile=$ARGV[1];
    open(NHANDLE,">$outfile") || die "error opning output file $outfile\n";
    #open(NHANDLE,">t");
    my $count=0;

    foreach (@NCLASS)
    {
       print STDOUT $mapSeq[$count]->[0]," ",$count," ",$mapSeq[$count]->[1],"\n";
       print NHANDLE "$_\n";
       $count++;
    }
    close NHANDLE;
       
0;          
