#!/usr/local/bin/perl

select STDERR; $| = 1; select STDOUT;

$minclustersize = 2;
$linkage = 1;

$usage .= "$0 - build clusters from MSPcrunch output\n";
$usage .= "\n";
$usage .= "Usage: $0 [-seqfile seqfile] [-verbose] [-minlen] [-minclustersize size] [-linkage linkage] <score cutoffs in order of increasing stringency>\n";
$usage .= "\n";
$usage .= "Program takes MSPcrunch input on STDIN and prints clusters to STDOUT\n";
$usage .= "\n";
$usage .= "Use \"/\" character to separate score, %ID, length cutoffs, e.g. \"100/60/.5\" means \"score >= 100, ID >= 60%, length >= .5 of longest sequence\"\n";
$usage .= "\n";
$usage .= "Use -seqfile option to specify a sequence file (required for lengths)\n";
$usage .= "Use -minlen to use length of shortest sequence, not longest\n";
$usage .= "\n";

while (@ARGV) {
    last unless $ARGV[0] =~ /^-/;
    $opt = lc shift;
    if ($opt eq "-seqfile") { defined($seqfile = shift) or die $usage; push @seqfile, $seqfile }
    elsif ($opt eq "-minlen") { $minlen = 1 }
    elsif ($opt eq "-minclustersize") { defined($minclustersize = shift) or die $usage }
    elsif ($opt eq "-linkage") { defined($linkage = shift) or die $usage }
    elsif ($opt eq "-verbose") { $verbose = 1 }
    else { die "$usage\nUnknown option: $opt" }
}

@ARGV>=1 or die $usage;
foreach (@ARGV) {
    ($score,$id,$lenfrac) = split /\//;
    push @score, $score;
    push @id, $id;
    push @lenfrac, $lenfrac;
    if ($lenfrac > 0) { $needlen = 1 }
}

foreach $seqfile (@seqfile) {
    warn "Reading sequence lengths from $seqfile...\n" if $verbose;
    open SEQFILE, $seqfile or die $!;
    $sep = $/;
    $/ = ">";
    $dummy = <SEQFILE>;
    while ($/ = "\n", $seqname = <SEQFILE>) {
	$seqname =~ s/\s.*//;
	$/ = ">";
	$seq = <SEQFILE>;
	chomp $seq;
	$seq =~ s/\s//g;
	$seqlen{$seqname} = length $seq;
    }
    $/ = $sep;
    close SEQFILE;
}

warn "Clustering...\n" if $verbose;

$sets = 0;
while (<STDIN>) {
    next if /^BLAST ERROR/ || !/\S/;

    @f = split;
    unless (defined($n1 = $index{$f[4]})) { push @name, $f[4]; $n1 = $index{$f[4]} = $#name }
    unless (defined($n2 = $index{$f[7]})) { push @name, $f[7]; $n2 = $index{$f[7]} = $#name }

    next if $n1 == $n2;

    unless ($seen{$f[4]}) {
	print STDERR "$f[4] " if $verbose;
	$seen{$f[4]} = 1;
    }

    if ($needlen) {
	exists($seqlen{$f[4]}) or die "Don't know length of $f[4]: $_";
	exists($seqlen{$f[7]}) or die "Don't know length of $f[7]: $_";
	if ($minlen) { $effseqlen = $seqlen{$f[4]} < $seqlen{$f[7]} ? $seqlen{$f[4]} : $seqlen{$f[7]} }
	else { $effseqlen = $seqlen{$f[4]} > $seqlen{$f[7]} ? $seqlen{$f[4]} : $seqlen{$f[7]} }
	$msplen = abs($f[3] - $f[2]) > abs($f[6] - $f[5]) ? abs($f[3] - $f[2]) : abs($f[6] - $f[5]);
    } else {
	$msplen = $effseqlen = 0;
    }

    for ($level=0;$level<@score;$level++) { last unless $f[0] >= $score[$level] && $f[1] >= $id[$level] && $msplen >= $effseqlen * $lenfrac[$level] }

    $oldlevel = $pairlevel{$pairkey = "$n1 $n2"};
    next unless $level > $oldlevel;

    delete $pairflag[$oldlevel]->{$pairkey} if $oldlevel;
    $pairflag[$level]->{$pairkey} = 1;
    $level{$pairkey} = $level;
}

print STDERR "\n" if $verbose;

@clustered = map(0,0..$#name);

for ($level=@score;$level>0;$level--) {
    %partners = ();
    foreach $pairkey (keys %{$pairflag[$level]}) {
	($n1,$n2) = split /\s/, $pairkey;
	next if $clustered[$n1] || $clustered[$n2];
	push @{$partners{$n1}}, $n2;
	push @{$partners{$n2}}, $n1;
    }

    foreach $n1 (reverse sort { @{$partners{$a}} <=> @{$partners{$b}} } keys %partners) {
	next if $clustered[$n1];
	$clustered[$n1] = 1;
	my @cluster = ($n1);
	for ($i=0;$i<$linkage;$i++) {
	    foreach $n (@cluster) {
		foreach $n2 (@{$partners{$n}}) {
		    push @cluster, $n2 unless $clustered[$n2];
		    $clustered[$n2] = 1;
		}
	    }
	}
	push @clusters, \@cluster if @cluster >= $minclustersize;
    }

}

foreach $cluster (sort { @$b <=> @$a } @clusters) { print join(" ",map($name[$_],@$cluster)),"\n" }
