#!usr/bin/perl

#########################################################
#
#Ian Odell
#Control Clustalw from Perl Code
#
#########################################################

use strict;
use warnings;
use Bio::DB::GenBank;
use Bio::Seq;
use Bio::SeqIO;
use Bio::SearchIO;
use Bio::Tools::Run::Alignment::Clustalw;
use Bio::AlignIO;
use Bio::Tools::Run::RemoteBlast;

####MAIN PROGRAM####
my $sequence = 'AY894554';
my $sequencefile = 'sequence.fasta';
my $db1 = 'nr';
my $clustal_dir = 'C:\Documents and Settings\Carl\Desktop\Perl\Clustalw_Perl';

my $seq_object  = dwnload($sequence, $sequencefile);
my $blast_filename = Blast_Seq($sequencefile, $db1);
#my $blast_filename = 'AY894554.blast';
my @seq_array = parseBlast($blast_filename);
runClustal(@seq_array, $sequence. "_aligned.aln", $clustal_dir);


###############################################################################################

#download a sequence file from a database of your choice
sub dwnload {
  my $sequence = $_[0];
  my $filename = $_[1];

  #create database object
  my $db = Bio::DB::GenBank->new();

  #get GenBank sequence
  my $seq_object = $db->get_Seq_by_acc($sequence); #get_Seq_by_id retrieves by GI number

  my $seqio_obj = Bio::SeqIO->new(-file => ">".$filename,
			          -format => 'fasta' );
  $seqio_obj->write_seq($seq_object);

  return ($seq_object);
}

###############################################################################################
  sub Blast_Seq {
    my $sequencefile = $_[0];
    my $db = $_[1];

    my $prog = 'blastn';
    my $e_val = '1e-10';

    my @params = (-prog => $prog,
		  -data => $db,
		  -expect => $e_val,
		  -readmethod => 'SearchIO');

    #$v is to turn on and off the messages
    my $v = 1;

    my $factory = Bio::Tools::Run::RemoteBlast->new(@params);

    my $str = Bio::SeqIO->new(-file => $sequencefile,
                              -format => 'Fasta');
    my $filename;
    while (my $input = $str->next_seq()) {
      my $r = $factory->submit_blast($input);
      print STDERR "Submitted Blast\nwaiting..." if ($v > 0);
      while (my @rids = $factory -> each_rid) {
	foreach my $rid (@rids) {
	  my $rc = $factory -> retrieve_blast($rid);
	  if (!ref($rc)) {
	    if ($rc < 0) {
	      $factory->remove_rid($rid);
	    }
	    print STDERR "." if ($v > 0);
	    sleep 5;
	  } else {
	    my $result = $rc-> next_result();
	    #save the output
	    $filename = $result -> query_name(). "\.blast";
	    $factory -> save_output($filename);
	    $factory -> remove_rid($rid);
	    print "\nQuery Name: ", $result -> query_name(), "\n";
	    while (my $hit = $result -> next_hit()) {
	      next unless ($v > 0);
	      print "\thit name is ", $hit->name, "\n";
	      while (my $hsp = $hit->next_hsp) {
		print "\t\tscore is ", $hsp->score, "\n";
	      }
	    }
	  }
	}
      }
    }
    return($filename);

  #non - OO way
  #uses the default database - nr in this case
  #my $blast_result = blast_sequence($seq);
  #write_blast($filename,$blast_result);
  }

###############################################################################################
sub parseBlast {
  my $filename = $_[0];

  # Get the report
  my $in = new Bio::SearchIO(-format => 'blast',
			     -file   => $filename);

  my @seq_array;
  my $count = 0;

  #parse through blast report and put all sequences into seq_array as seq objects
  while ( my $result = $in->next_result ) {
    while ( my $hit = $result->next_hit ) {
      while ( my $hsp = $hit->next_hsp ) {
	push(@seq_array, Bio::Seq->new(-seq => $hsp->hit_string,
				       -id => $hit->accession. "_$count"));
	$count++;		#make sure unique seq id

      }
    }
  }

  #print the id and sequence of each seq object in @seq_array
  #and write all sequences to a fasta file
  my $i;
  my $temp_str;
  my $seq_io = Bio::SeqIO->new(-file => '>blast2fasta.fasta', #>> means append to file
			       -format => 'Fasta');
  for $i (@seq_array) {
    $seq_io->write_seq($i);
    $temp_str.=">". $i->id. "\n". $i->seq. "\n\n";
  }
  print $temp_str;
  return (@seq_array);
}


################################################################################################

#Run clustal from this program (I think Clustalw has to be in same directory as perl file)
sub runClustal{

  my $prog_dir = pop(@_);
  my $filename = pop(@_);
  my @seq_array = @_;

  #align @seq_array with clustalw
  my @params = ('ktuple' => 3,
		'matrix' => 'BLOSUM');

  my $factory = Bio::Tools::Run::Alignment::Clustalw->new(@params);

  $factory->program_dir($prog_dir);

  #make reference to an array of Bio::Seq objects
  my $seq_array_ref = \@seq_array;

  my $aln = $factory->align($seq_array_ref);

  #write alignment to file

  my $out = Bio::AlignIO->new(-file => ">". $filename,
			      -format => 'clustalw');
  $out->write_aln($aln);
}
