#!/usr/bin/perl -w
# ANALYZE_SNPS 
#    VERSION: Version 1 (29 June 2004)
#    PURPOSE: Accumulate data on how close SNPs are to CG dinucleotides
#    INPUT FILE:  
#             Sequence
#             Sequence
#             . . .
#             Sequence
#
#       Each sequence consists of
#             5' region (up to 300 nucleotides)
#             Variable nucleotide (R Y S W K M)
#             3' region (up to 300 nucleotides)
#
#       Each sequence lies on a single long line
#
#    OUTPUT FILE: 
#             

############## LIBRARIES AND PRAGMAS ################

  use strict;

#################### CONSTANTS ######################

  my $true = 1;                # Perl has no specific value for true
                               #     so we supply one
  my $false = not $true;       # Equivalent to $false = 0
  my $LF = "\n";               # Linefeed
  my $tab = "\t";              # Tab
  my $max_length = 300;

#################### VARIABLES ######################

  my $sequence;                   # Current sequence to be analyzed
  my $sequences_analyzed = 0;     # Number of sequences processed
  my $total_nucleotides = 0;      # Number of nucleotides processed
  my @count = (0) x $max_length;  # Number of SNPS at position N away from CG
  my $left_length;                # Length from CG to variable position 5'->3'
  my $right_length;               # Length from variable position to CG 5'->3'
  my $length_to_CG;               # Minimum of $left_length and $right_length

###################### FILES ########################

  my $input_file = "rs_ch21.txt";
  open INPUT_FILE, "<$input_file" 
                        or die "Can't open $input_file: $!\n";

  my $output_file = "rs_vs_cg.txt";
  open OUTPUT_FILE, ">$output_file" 
                        or die "Can't open $output_file: $!\n";

################### MAIN PROGRAM ####################

  $sequence = <INPUT_FILE>;

  while ($sequence) {
     chomp $sequence;
   
     if ($sequence =~ ### Put in regular expression here ###)
        {$left_length = length ($1)}
     else 
        {$left_length = $max_length}
     if ($sequence =~ ### Put in regular expression here ###)
        {$right_length = length ($1)}
     else 
        {$right_length = $max_length}

     $length_to_CG = Min ($left_length, $right_length);
     if ($length_to_CG < $max_length) {
        ### Put something here to increment array count by 
             one at position of $length_to_CG.
             It will start either:
                $count...      = $count...     + 1;
             or
                ++$count...
        }

     $sequences_analyzed = $sequences_analyzed + 1;
     $total_nucleotides = $total_nucleotides + length($sequence);
     $sequence = <INPUT_FILE>;
  }     

  print $sequences_analyzed, " sequences analyzed", $LF;
  print $total_nucleotides, " nucleotides analyzed", $LF, $LF;
  Save_counts();

  print "Contents of array:", $LF;
  print join (" ", @count), $LF;

  print "ALL DONE!";

#################### SUBROUTINES ####################

#### SAVE_COUNTS
#
  sub Save_counts {
    foreach $length_to_CG (1 .. $max_length-1) {
       print OUTPUT_FILE $length_to_CG, $tab, $count[$length_to_CG], $LF;
    }
  }

#### MIN (m and n)
#    Returns the lesser of the two numbers
  sub Min {
     my ($a, $b) = @_;
     if ($a < $b) { return $a }
     else { return $b }
  }
