#!/usr/bin/perl -w
use strict;

#############################################
# Authors: Marta Massanella & Ingrid Pallās #
# Date: March 2005                          #
#############################################



### This program needs the output file from findsnp.pl ###

if ( scalar(@ARGV) < 2 ) {
  print STDERR "secis.pl <SequenceFile> <PolymorphismFile>\n";
  exit(1);
}

open ( MYSEQ, "<$ARGV[0]" );                                 # Secis Sequence from the Selenporotein analyzed in findsnp.pl is required
my @seq = <MYSEQ>;
close ( MYSEQ );


open ( MYSNP, "<$ARGV[1]" );                                 # Output file from findsnp.pl is required.
my @SNP = <MYSNP>;
close ( MYSNP );

### MY VALUES ###
my $i;
my $j;

my @seq2 = ();


##############################################################################
#     PROGRAM BODY: Introducing one SNP in a Secis query sequence            #
##############################################################################


@seq2 = split ( //,$seq[1] );                               # The vector where is the input file (@seq) has 2 positions:it is splitted in
                                                            # in @seq2 the nucleotides from the  seq which it is found in the second
                                                            # position of the vector ($seq[1]).
open( FILE, ">seqSecis.txt" );                              # All the seq that we will obtain we will be stored in seqSecis.txt.

$j=0;                                            	    # $j is an index which allows us to scan the whole input using an iterative composition
while ( $j < scalar(@SNP) ) {

    if ( $SNP[$j] =~ /^SNP\s+\w+\:\s+(\d+)([acgt\-])(.*?)\s+;/o ) {  # It is select only the SNPs from the input file & we store the
                                                                     # position ($position),the nucleotid of the query sequence ($query)
	my $position = $1 - 1;                                       # & the change/s from subject sequence (@change). This last vector
        my $query = $2;                                              # is splitted in order to have each nucleotide change in each
	my @change = split ( /\//, $3 );                             # position of the vector.

	    my $k = 0;                                            # $k is an index which allows us to scan the whole changes from one seq
	    while ( $k < scalar(@change) ) {


		$i=0;                                             # $i is an index which allows us to scan each nucleotide of secis seq
		while ( $i <scalar(@seq2) ) {

		    if ( $i == $position ) {	                  # If the position from the SNP coincide with $i which is scanning the
			                                          # seq, two different situations can occur:
			if ( $query eq '-' ) {                    # 1. We need to introduce a nucleotide to Query seq because it has
			    my $temp = $change[$k];               # successed a deletion in query seq or an insertion in the EST.
			    
				print FILE ">SECISseq",$position+1,",-/$temp\n".join('',@seq2[(0..$i)],$temp,@seq2[($i+1..$#seq2)])."\n";
			}                                         
								  # The change is kept in a temporal vector, & afterwards it is printed
                                                                  # the seq until $i position, the change and finally the following seq
                                                                  # ($i+1 until the end).
			else {
			    $change[$k] = uc($change[$k]);        # 2.We need to change a nucleotide from query seq at position $i.First
			    my $seq3 = $seq2[$i];                 # at all, it is kept the nucleotide from query seq,afterwards we change
			    $seq2[$i] = $change[$k];              # this nt to the corresponding nt & it is print the result seq.
			    
			    print FILE ">SECISseq", $position+1,",$query/$change[$k]\n".join('',@seq2)."\n";
			    $seq2[$i] = $seq3;                   
								  # The old nt is replaced again in query seq after this step.
			}
		    }


		    $i = $i +1;
		}


		$k = $k + 1;
	    }
    }

    $j = $j + 1;
}

close (FILE);


print "Result in seqSecis.txt file\n";