#!usr/bin/perl
######################################################################
#
#Ian Odell
#Read in alignment file
#Find the beginning and end of each sequence
#Fill in any dashes with a standard sequence for each type
#
######################################################################

use strict;
use warnings;
use Bio::Seq;
use Bio::SeqIO;
use Bio::PrimarySeq;

my $inputfilename = "HLA_ABC_fix_rem.fasta";
my $outputfilename = ">output.fasta";

#read in input fasta file
my $in = Bio::SeqIO->new(-file => $inputfilename);
my @input_seq_array;
while (my $seq = $in->next_seq()) {
  push @input_seq_array, $seq;
}

#define a default sequence with which the other sequences will fill in
my $gapA = $input_seq_array[0]->seq;
my $gapB = $input_seq_array[42]->seq;
my $gapC = $input_seq_array[97]->seq;

# for each sequence, find the first and last base
# replace the dashes at the beginning and end with the profile sequence
my $sequence;
my $begin;
my $end;
my $gap;
my $name;
my $output;
foreach my $seq (@input_seq_array){
  $sequence = $seq->seq;
  $begin = BaseIndex($sequence);  #find first base
  $end = RBaseIndex($sequence);  #find last base
  $name = $seq->id;
  $gap = "";
  #change gap according to type
  if ($name =~ 'A'){
  $gap = $gapA;
  } elsif (substr($name, 0, 1) eq 'B'){
  $gap = $gapB;
  } elsif (substr($name, 0, 1) eq 'C'){
  $gap = $gapC;
  } else {}
  substr($sequence, 0, $begin, substr($gap, 0, $begin));  #replace dashes with gap string
  substr($sequence, $end+1, length($sequence)-($end+1), substr($gap, $end+1));  #replace end dashes
  $output .= ">" . $seq->id . "\n$sequence\n";
}
#open output file and write PWM to it
open FH, $outputfilename or die $!;
print FH $output;

###################subroutine to find the first base########################
sub BaseIndex{
  my $string = $_[0];
  my $baseindex = -1;
  my $index = 0;
  my $bases = 'ACGT';
  while ($baseindex == -1 && $index < length($string)) {
    if ($bases =~ substr($string, $index, 1)) {
      $baseindex = $index;
    } else {
      $index++;
    }
  }
  return $baseindex;
}
####################subroutine to find the last base#########################
sub RBaseIndex{
  my $string = $_[0];
  my $baseindex = -1;
  my $index = length($string)-1;
  my $bases = 'ACGT';
  while ($baseindex == -1 && $index >= 0) {
    if ($bases =~ substr($string, $index, 1)) {
      $baseindex = $index;
    } else {
      $index--;
    }
  }
  return $baseindex;
}
#############################################################################
