#! -w
#
# pick-key-residues  
#	[-minbits 2.0]	[-minfreq 0.01]
#	 < x.saves >  foo.res
#
# selects those residues in the .saves file 
#	that have bit-savings (relative entropy) of minbits or more
#	and outputs them in a format suitable for use by select-by-key-residues
# The amino acids output for each key residue are those with at least
#	minfreq probability on the appropriate line of the .saves file.
#
# INPUT format:
#	A .saves file has one line per alignment column of a multiple
#	sequence alignment.
#	Each line has 41 fields:
#		bit savings
#		20 (amino acid, probability)
#	with the amino acids sorted into order by decreasing probability.
#
# OUTPUT format:
# The format for the foo.res file is a white-space separated
# list of residues that must be present.
# For example,
#	N150 
#	AS180 
#	P-225
# would mean that alignment column 150 (with 1 as the first column)
# must have Asparagine, 180 must have Alanine or Serine, and 225 must
# have Proline or a gap.
# (Gaps are not output by this program.)

{   use FileHandle;
    STDOUT->autoflush(1);
    STDERR->autoflush(1);
    use English;
    use Getopt::Long;

    $min_bits = 2.0;
    $min_freq= 0.1;
    GetOptions("minbits=f" => \$min_bits,
    	"minfreq=f" => \$min_freq);
	
    $column=0;
    while (<STDIN>)
    {   s/^\s*//;	# remove leading blanks
        my ($saves, @fields) = split(/\s/);	# split into fields
    	next if (scalar(@fields)==0);
	$column++;
	if ($saves >= $min_bits)
	{   my $aastring="";
	    for (my $i=0; $fields[2*$i+1]>=$min_freq; $i++)
	    {    $aastring .= $fields[2*$i];
	    }
	    if (length($aastring)==0)
	    {    print STDERR "Warning: column $column has $saves bits saved,"
	    		. " but no amino acid with frequency >= $min_freq\n";
	    }
	    else
	    {	print "$aastring$column\n";
	    }
	}
    }
}
