#!/usr/bin/perl -w

# script to construct single-linkage clusters
#
# input:  file with two columns, each entry one sequence id, each line indicating a "hit",
#	Format of the sequence ids is arbitrary, as long as each entry is followed by white space.
#	This script can be used for any input where each line is an edge and the pair of words are nodes.
#	
# output: file with two columns, first column cluster number, in "clXX" format, second column sequence id.
#	in the same format as in the input file.
#
# strategy:  build data structures "on the fly" (i.e., as the data is read into memory),
#	linking clusters as necessary, making sure to merge smaller into larger to save time.
#	
# Reference:  This is based on the algorithm developed in "blink.c" by MJ Sanderson, which see for 
#	further documentation.  More recently, MJS has produced "blinkcomp.c" which employs a 
#	linear-time algorithm to perform the clustering using connected components.  All three 
#	programs take the same input format and produce the same basic output, but have
# 	slight differences in the way the data is summarized, e.g., in the summary file here.
#
$Usage = "

Usage:
		./blinkPerl.pl -i infile -o outfile -s summaryFile [-v logFile]

";
#	example:
#		./blinkPerl.pl -i dna5k.hitlist -o dna5k.clList -s dna5k.clSum -v dna5k.blinkLog
#
$version = "\nblinkPerl.pl v. 1.0
  M. McMahon 14 Mar 05\n";
print "$version";
#
# command line input ---------------------------------------------------
# ----------------------------------------------------------------------

while ($fl = shift @ARGV)
  {
  $par = shift @ARGV;
  die ("$Usage") if (! $par);
  if ($fl =~ /-i/) {$infile = $par;}
  if ($fl =~ /-o/) {$outfile = $par;}
  if ($fl =~ /-s/) {$sumfile = $par;}
  if ($fl =~ /-v/) {$logfile = $par;}
  }
die ("$Usage") if (! ($infile && $outfile && $sumfile));
die ("Input file not available.\n") if (! -e $infile);
die ("Output file already exists.\n") if (-e $outfile);
die ("Summary file already exists.\n") if (-e $sumfile);
if ($logfile) {$verbose=1;}
else {$verbose=0;}

# initialize data structures ----------------------------

%siClH = ();    # keys are seqids pointing to cluster numbers
%clSiHoA = ();  # keys are cluster numbers pointing to arrays of seqids

# -----------------------------

$nextCl = 1;

open FH1, "<$infile";
if ($verbose) {open FH3, ">$logfile";}
while (<FH1>) 
  {
  chomp;
  ($si1, $si2) = /\S+/g; 
  # if the two are equal, start a new cluster if necessary
  if ($si1 eq $si2)  
    { 
      if (! (exists $siClH{$si1})) 
        {
        $siClH{$si1}=$nextCl; 
        $clSiHoA{$nextCl} = [ $si1 ]; 
        ++$nextCl;
        next;
        }
    }
  # if both seqids have been encountered, the two clusters merge, smallest merging into largest
  if (exists $siClH{$si1} && exists $siClH{$si2}) 	
    {
    ($cl1, $cl2) = ($siClH{$si1}, $siClH{$si2});
    next if ($cl1 == $cl2);
    $size1 = scalar @{$clSiHoA{$cl1}};
    $size2 = scalar @{$clSiHoA{$cl2}};
    if ($size1 < $size2) {$oldCl = $cl2; $newCl = $cl1;}
    else {$oldCl = $cl1; $newCl = $cl2;}
    if ($verbose) {print FH3 "merging cluster $oldCl into $newCl because of $si1 or $si2\n";}
    for $si (@{$clSiHoA{$oldCl}}) 
      {$siClH{$si} = $newCl;}
    push @{$clSiHoA{$newCl}}, @{$clSiHoA{$oldCl}};
    delete $clSiHoA{$oldCl};
    }
  # if only one has been encountered, add the other seqid to the cluster 
  elsif (exists $siClH{$si1}) 
    {
    $cl = $siClH{$si1};
    $siClH{$si2} = $cl;
    push @{$clSiHoA{$cl}}, $si2;
    } 
  elsif (exists $siClH{$si2}) 
    {
    $cl = $siClH{$si2};
    $siClH{$si1} = $cl;
    push @{$clSiHoA{$cl}}, $si1;
    }
  # if neither is listed, start a new cluster.
  else
    {
    if ($verbose) {print FH3 "making new cluster $nextCl of $si1 and $si2\n";}
    $siClH{$si1} = $nextCl;	
    $siClH{$si2} = $nextCl;
    $clSiHoA{$nextCl} = [ $si1, $si2 ];
    ++$nextCl;	 	
    }
  }
close FH1;

# remove dups from cluster lists ------------------

%globalSiH = ();  # this is to make sure that no seqid is in 2 dif clusters
$globalSiHref = \%globalSiH;

for $cl (keys %clSiHoA)
  {
  removeDups ($clSiHoA{$cl}, $globalSiHref);
  $length2 = scalar @{$clSiHoA{$cl}};
  die ("cl $cl has no seqids\n") if ($length2 == 0);
  }

# rename clusters with sequential numbers ---------
#	(creating new clHoA, old structures remain the same)

$newCl = 1;
for $oldCl (keys %clSiHoA)
  {
  $newClSiHoA{$newCl} = $clSiHoA{$oldCl};
  ++$newCl;
  }

# print output ------------------------------------
@allSis = (keys %siClH);
@allCls = sort {$a <=> $b} (keys %newClSiHoA);
$totalCls = scalar @allCls;
$totalSis = scalar @allSis;
$checkTot = 0;

open FH4, ">$sumfile";
print FH4 "Total number of sequences: $totalSis.\nTotal number of clusters: $totalCls.\n";
print FH4 "\nCluster\tsequences\n";
open FH2, ">$outfile";

for $cl (@allCls)
  {
  for $si (@{$newClSiHoA{$cl}})
    {
    print FH2 "cl$cl\t$si\n"; 
    }
  $numSis = scalar @{$newClSiHoA{$cl}};
  die ("NO seqids in cl $cl") if ($numSis == 0);
  print FH4 "$cl\t$numSis\n";
  $checkTot = $checkTot + $numSis;
  }
  {print FH4 "Number of seqids assigned to clusters is $checkTot and total seqids in was $totalSis\n";}
close FH2;
close FH3;
close FH4;

# -------------------------------------------------
# sub removeDups ----------------------------------
#  takes ref to an array
#	modifies original array
#	also takes ref to global si hash and adds to it...
#  returns nothing

sub removeDups
{
my ($aref, $href) = @_;
my @oldA = @{$aref};
my %tempH = ();
my @newA = ();

for $x (@oldA)
  {
  if (! (exists $tempH{$x}))
    {
    if (! (exists ${$href}{$x}))
      {${$href}{$x}=1;}
    else {die ("$x is in mulitple clusters\n");} 
    push @newA, $x; 
    $tempH{$x} = 1;
    } 
  }
@{$aref} = @newA;
return;
}

