#!usr/bin/perl
######################################################################
#
#Ian Odell
#Read in alignment file of unique sequences - HCV-Study#_type_duplicates
#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::String;


my $inputfilename = "type6_ambigsRemoved.aln";
my $outputfilename = ">type6PWM_ambigsremoved.txt";

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

#calculate total number of sequences of this type
my $total=0;
foreach my $seq ($aln->each_seq){
  $total += substr($seq->id, rindex($seq->id, '_')+1);
}
print $total, "\n";

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 $numA;
my $numC;
my $numG;
my $numT;
my $numDash;
my $numAmbig;
my $format = new Number::Format();
my $id;

while ($index <= $aln->length) {
  $numA = 0;
  $numC = 0;
  $numG = 0;
  $numT = 0;
  $numDash = 0;
  $numAmbig = 0;
  foreach my $seq ($aln->each_seq) {
    $res = $seq->subseq($index, $index);
    $id = $seq->id;
    #if the residue is A, C, G, T, -, or Ambiguous, add the number of duplicates
    if ($res eq 'A') {
      $numA += substr($id, rindex($id, '_')+1);
    } elsif ($res eq 'C') {
      $numC += substr($id, rindex($id, '_')+1);
    } elsif ($res eq 'G') {
      $numG += substr($id, rindex($id, '_')+1);
    } elsif ($res eq 'T') {
      $numT += substr($id, rindex($id, '_')+1);
    } elsif ($res eq '-') {
      $numDash += substr($id, rindex($id, '_')+1);
    } else {			#ambiguous bases: R, Y, K, M, S, W, B, D, H, V, N
      $numAmbig += substr($id, rindex($id, '_')+1);
    }
  }
#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, 7);
  push @pC, $format->round($numC/$total, 7);
  push @pG, $format->round($numG/$total, 7);
  push @pT, $format->round($numT/$total, 7);
  push @Dash, $format->round($numDash/$total, 7);
  push @Ambig, $format->round($numAmbig/$total, 7);
  $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";
}

my $all = "$ASTRING\n" . "$CSTRING\n" . "$GSTRING\n" . "$TSTRING\n" . "$DashSTRING\n" . "$AmbigSTRING\n";

open FH, $outputfilename or die $!;
print FH $all;


