#!/usr/bin/perl
use strict;
use warnings;

sub iterblast {
    # settings
    my $myblast = "blastall -p blastn -m 8";
    my $starterseq =
        "/share/home/tjuh/Thema12-2011/Project/testblast/matchedhits.fasta";
    my $blastdb =
        "/share/home/tjuh/Thema12-2011/Project/databases/termite_hindgut/hindgut_db";
    my $outputroot = "/share/home/tjuh/Thema12-2011/Project/testblast/test1";
    my $evalue     = 0.01;
    my $wordsize   = 11;
    my $roundlimit = 4;

    # output files:
    my $seq_db_fasta = $outputroot . ".tmp";
    my $seq_db_all   = $outputroot . ".all";
    my $logfile      = $outputroot . ".log";

    # remove existing files
    unlink $seq_db_fasta;
    unlink $seq_db_all;
    unlink $logfile;

    open my $log, '>', $logfile or die "Can't open $logfile: $!\n";
    print $log <<LOGMSG;
Configurations:
Input file       : $starterseq
Output files     : $outputroot\-
Blast Engine     : $myblast
Target database  : $blastdb
E-value cut-off  : $evalue
Blast word size  : $wordsize
Blast round limit: $roundlimit

Runs:
LOGMSG

    # hits hash
    my %hits;
    my %queryids;

    # Start the iterative blastn rounds
    for my $round (1 .. $roundlimit) {
        my $saveblst = $outputroot . ".txt";
        
        unlink $saveblst;

        # print $log info
        my $curtime = localtime();
        my $infile = $round == 1 ? $starterseq : $seq_db_fasta;

        print $log "$curtime - Blast seq: $infile\n";
        print $log "$curtime - Blast db : $blastdb\n";

        # run blast
        system(
            "$myblast -i $infile -d $blastdb -o $saveblst -e $evalue -W $wordsize"
            );

        open my $parse, '<', $saveblst or die "Can't open $saveblst: $!\n";
        open my $fasta, '>', $seq_db_fasta
            or die "Can't create $seq_db_fasta: $!\n";

        # parse blast output (tab form) and add matched id's in hits hash
        while (defined (my $line = <$parse>)) {
            chomp $line;
            my @fields = split(/\t/, $line);
            ++$queryids{$fields[0]};
            ++$hits{$fields[1]};
        }

        # get header and sequence with 'fastacmd' and write to new fasta file
        my $hitcount;

        while (my $key = each %hits) {
            next if exists $queryids{$key};

            # only write uniq ids to file
            system("fastacmd -d $blastdb -s 'lcl|$key' >> $seq_db_fasta");
            $hitcount++;
        }

        close $fasta;
        close $parse;

        if ($hitcount) {
            print $log "Round $round has $hitcount unique hits\n\n";
            next;
        }

        print $log "No hits in Round $round ==> STOP\n";
        last;
    }

    print "Total uniq hits found in all rounds: ", scalar keys %hits, "\n";
    close($log);
}