#!usr/bin/perl
######################################################################
#
#Ian Odell
#Read in alignment file
#Output Position Weight Matrix
#
######################################################################

use strict;
use warnings;

use Bio::AlignIO;
use Bio::Align::AlignI;
use Bio::Matrix::PSM::PSM;
use Bio::Matrix::PSM::SiteMatrix;
use Number::Format;

#use IO::File::String;

my $inputfilename = "HLA_ABC_exon3.aln";
my $outputfilename = ">HLA_ABC_exon3_PWM.txt";

my $in  = Bio::AlignIO->new(-file => $inputfilename);

my $aln = $in->next_aln();

my @pA;
my @pC;
my @pG;
my @pT;
my @Dash;
my @Ambig;

#For each index in the alignment, calculate the frequency of each base
#and add to frequency arrays
my $index = 1;
my $res;
my $total = $aln->no_sequences;
my $numA;
my $numC;
my $numG;
my $numT;
my $numDash;
my $numAmbig;
my $format = new Number::Format();

my $seq;
while ($index <= $aln->length) {
  $numA = 0;
  $numC = 0;
  $numG = 0;
  $numT = 0;
  $numDash = 0;
  $numAmbig = 0;
  foreach $seq ($aln->each_seq) {
    $res = $seq->subseq($index, $index);
    if ($res eq 'A') {
      $numA++;
    } elsif ($res eq 'C') {
      $numC++;
    } elsif ($res eq 'G') {
      $numG++;
    } elsif ($res eq 'T') {
      $numT++;
    } elsif ($res eq '-') {
      $numDash++;
    } else {			#ambiguous bases: R, Y, K, M, S, W, B, D, H, V, N
      $numAmbig++;
    }
  }
#before pushing the values onto their respective array, format value and divide
#by total number of sequences for percentage
  push @pA, $format->round($numA/$total, 6);
  push @pC, $format->round($numC/$total, 6);
  push @pG, $format->round($numG/$total, 6);
  push @pT, $format->round($numT/$total, 6);
  push @Dash, $format->round($numDash/$total, 6);
  push @Ambig, $format->round($numAmbig/$total, 6);
  $index++;
}

my $n;
my $ASTRING = "A\t";
foreach $n (@pA){
  $ASTRING = "$ASTRING" . "$n\t";
}
my $CSTRING = "C\t";
foreach $n (@pC){
  $CSTRING = "$CSTRING" . "$n\t";
}
my $GSTRING = "G\t";
foreach $n (@pG){
  $GSTRING = "$GSTRING" . "$n\t";
}
my $TSTRING = "T\t";
foreach $n (@pT){
  $TSTRING = "$TSTRING" . "$n\t";
}
my $DashSTRING = "DASH\t";
foreach $n (@Dash){
  $DashSTRING = "$DashSTRING" . "$n\t";
}
my $AmbigSTRING = "AMBIG\t";
foreach $n (@Ambig){
  $AmbigSTRING = "$AmbigSTRING" . "$n\t";
}

#concatenate the strings
my $all = "$ASTRING\n" . "$CSTRING\n" . "$GSTRING\n" . "$TSTRING\n" . "$DashSTRING\n" . "$AmbigSTRING\n";

#open output file and write PWM to it
open FH, $outputfilename or die $!;
print FH $all;
