#!/usr/local/bin/perl -w
use strict;

############################# Global variables
#
my @sequences;           # The set of sequences in which we search for
                         #    a consensus sequence.

my $sequence_length;     # The length of each sequence. (We make the slightly
                         # artificial restriction that all sequences be of the
                         # same length.)

my @consensus;           # The consensus sequence.
my $highThreshold = 0.8; # Print consensus base in upper case if this fraction
                         #    of sequences agree.
my $lowThreshold = 0.5;  # Print consensus base in lower case if this fraction
                         #    of sequences agree.  (Must be at least 0.5, and
                         #    less than $highThreshold.)

my @bases = ("G", "A", "T", "C");

############################ Main program
#
get_sequences('NtcA-binding.txt');
calculate_consensus_sequence();
print_sequences();
print_consensus();

########################### Subroutines

#### get_sequences
#
# Call as
#
#      get_sequences($name_of_file)
#
# Fills the @sequences array

sub get_sequences {
   my ($sequencePath) = @_;
   open INPUT, "< $sequencePath" or die "Can't open $sequencePath: $!";
   while(<INPUT>) {
      next if /^#/;
      next if /^\s*$/;
      chomp;
      my ($sequence) = /^.*\t(.*)$/
         or die "Don't understand line $. of $sequencePath";
      push @sequences, [split("", uc($sequence))];
      $sequence_length = length($sequence);
   }
   if (@sequences == 0) {
       die "No sequences found in $sequencePath\n";
   }
   return @sequences;
}

#### print_sequences
#
# Print @sequences, one sequence per line
#
sub print_sequences {
   for (my $row = 0; $row < @sequences; $row = $row+1) {
      for (my $column = 0; $column < $sequence_length; $column = $column+1) {
         print $sequences[$row][$column];
      }
      print "\n";
   }
}

#### calculate_consensus_sequence
#
# For each column, set
#
#   $consensus[$column] to "A"  if at least $highThreshold fraction of bases
#                               in the column are "A"
#   $consensus[$column] to "a"  if the fraction of "A"'s in the column is
#                               at least $lowThreshold and less than
#                               $highThreshold.
# ... and similarly for T, G, and C.
#
# Set $consensus[$column] to "." if none of the bases meet those conditions.
#
sub calculate_consensus_sequence {
   for (my $column = 0; $column < $sequence_length; $column = $column+1) {
      my @count = (0, 0, 0, 0);
      for (my $row = 0; $row < @sequences; $row = $row+1) {
         for (my $b = 0; $b < @bases; $b = $b+1) {
            if ($sequences[$row][$column] eq $bases[$b]) {
               $count[$b] = $count[$b] + 1;
            }
         }
      }
      $consensus[$column] = ".";
      for (my $b = 0; $b < @bases; $b = $b+1) {
         if ($count[$b]/@sequences >= $highThreshold)
            { $consensus[$column] = $bases[$b] }
         elsif ($count[$b]/@sequences > $lowThreshold)
            { $consensus[$column] = lc($bases[$b]) }
      }
   }
}

#### print_consensus
# 
sub print_consensus {
   for (my $column = 0; $column < $sequence_length; $column = $column+1) {
      print $consensus[$column];
   }
   print "\n";
}
