#!/usr/bin/perl -w
#
#       SAM: Sequence Alignment and Modeling Software System                 
#  (c) 1992-1999 Regents of the University of California, Santa Cruz
#
#          Sequence Alignment and Modeling Software System
#         http://www.cse.ucsc.edu/research/compbio/sam.html
#
# ------------- Citations (SAM, SAM-T99, HMMs) -----------------
# R. Hughey, A. Krogh, Hidden Markov models for sequence analysis:
#  Extension and analysis of the basic method, CABIOS 12:95-107, 1996.
# K. Karplus, C. Barrett, R. Hughey, Hidden Markov models for detecting
#  remote protein homologies, Bioinformatics 14(10):846-856, 1999.
# A. Krogh et al., Hidden Markov models in computational biology:
#  Applications to protein modeling, JMB 235:1501-1531, Feb 1994.
# --------------------------------------------------------------
#
#          Distributed with signed license agreement only.
#  
#          Questions or comments to sam-info@cse.ucsc.edu
#


#!/usr/bin/perl -w

# target2k
# builds an alignment of probable homologs of an initial seed sequence
# or alignment, with homologs chosen from NR.
#
# This script is an updated and re-parameterized version of the older
# target99, target98-with-param, and target98 scripts.
#
# WARNING: the defaults may change.
#
# This script takes one required and several optional arguments:
#  REQUIRED:
#	-out root		creates root.a2m for output
#				(also root.cst if constraints are
# 				turned on)
#  OPTIONAL MAJOR CHOICES:
#	-seed file.a2m		a guide sequence (or alignment)
#				The seed is usually expected, but if none
#				is provided and -homologs or -close is
#				specified, then an initial alignment will
#				be built from the (close) homologs.
#	-superfamily		Sets many options appropriately for
#				modeling a superfamily given a single sequence.
#				[Default choice.  Tuned somewhat for SCOP test]
#	-family			Sets options appropriately for
# 				modeling a family.  These options have NOT
#				been carefully tuned, and are just a
# 				guess at initial parameters.
#	-tuneup			Eliminates search, assumes that all
# 				that is wanted is improvement of
# 				initial alignment (weeding out
# 				too-distant sequences, realigning).
#

# OPTIONAL MINOR CHOICES:
#	-iter 4			number of iterations
#	-all			creates root_1.a2m, root_2.a2m, ...
#					for the multiple iterations
#	-thresholds 0.001,0.02,0.2,1.0 specifies the E-value
#				thresholds to use for each iteration
#				(E-value is expected number of false positives)
#	-db_size  408000	Specify number of sequences in database to assume for
#				E-value calculations.  If not specified, then
#				size of database in -db is assumed.
#       -homologs               File containing potential homologs of the guide
#                               sequence. If this parameter is not specified
#                               the default action is to extract homologs
#                               from NR using blast or Smith-Waterman search.
#       -close                  File containing close homologs of the guide
#                               sequence. These close homologs are used to
#                               train the model on the first training
#                               iteration. If this parameter is not specified,
#                               but -homologs is specified, then the file from
#                               the -homologs option is used for the close
#                               homologs. If neither -close nor -homologs are
#                               specified, then the default is obtain both
#                               from NR using blast, but -close uses a
#				much tighter threshold.
#	-no_search		use only the -seed and -homologs or -close
#				homologs--don't look for more.
#				-db_size should probably be specified
#				with -no_search
#       -constraints 1          Should constraints be used?
#				0 means no constraints.
#				1 means use constraints on seed
# 					alignment only.
#				(Right now, constraints are for the entire
#				alignment, but this should be changed to
#				allow user-specified initial constraints.)
#	-force_seed 1		set (default) to force all of seed
#				into the training set.
#	-db nr			Which non-redundant database to use.
#				Note: database MUST be in FASTA format
#				and have a wu-blastp index.
#	-nsurgery 0		number of surgery steps to allow on
#				each iteration	(default 0)
#				(Note: surgery always allowed in
#				 creating initial alignment when no
#				 seed is provided.) 
#				-nsurgery >0 and -force_seed 1 are incompatible.
#	-no_sort		Turn off sorting of sequences.
#
#       -ncbiblast              (Default): Use NCBI blast prefilter
#       -wublast                Use WU-Blast prefilter
#

# VERY MINOR OPTIONS (developer use only?)
#	-tmp_dir /projects/compbio/tmp	where to build the temporary 
#					subdirectories
#	-jump_in  0.2		probability associated with jump-in and
#	-jump_out 1.0			jump-out on local alignment
#	-fimstrength 1.0	probability multiplier for letters in fims
#	-fimtrans 1.0		another parameter affecting loops on fims
#	-cut_match 0.2		In surgery, if fewer sequences than
#				this fraction use a match state, delete it.
#	-cut_insert 0.3		In surgery, if more sequences than this
#				fraction use an insert state, add match states.
#	-frac_insert 0.3	In surgery, replace insert state with
#				frac_insert*avg number of characters generated
#	-reverse_diff 4		how much tighter is the simple
# 				threshold than the reversed-sequence threshold
#				(This is just a heuristic to reduce
#			the number of reverse-scores that are computed.)
#	-prefilter_thresholds 1.e-3,20,400,3000
#				E_value thresholds used for each iteration
#	-blast_max_report 20000	Maximum number sequences reported by WU-BLAST
#
# 	-aweight_method 2	Henikoff weighting
#				(1=EntropyWeight, 2=HenikoffWeight, 3=FlatWeight)
#	-aweight_bits 0.5	(can say "var50" to get variable weighting
#				 50/predicted_alignment_length---can use
#				 comma-separated list to get different weight
#				 on each iteration)
#	-aweight_exponent 0.5	Remember to change to 10 or so if using
#				EntropyWeight (-aweight_method=1)
#	-mixture recode3.20comp Dirichlet mixture to use for
#				regularizing amino acid distributions.
#	-init_adpstyle 1	Alignment style to use for
#				building initial alignment from unaligned seqs
#	-search_adpstyle 1	Alignment style to use for selecting
#				training set for retraining HMM
#	-realign_adpstyle 1	Alignment style to use for making alignment
#				from trained HMM (internal)
#	-final_adpstyle 5	Alignment style to use for making alignment
#				from trained HMM (final)
# 	-init_trans_reg bl29.regularizer	transition regularizer
#				for first iteration
#	-mid_trans_reg long_match.regularizer	transition regularizer for
#				middle iterations
#	-final_trans_reg fssp-trained.regularizer transition regularizer for
#				final iteration
#	-keep_temporary 0	don't delete temporary directories
#	-full_seq_align 0	Use full sequence alignment, not
#				multiple-domain for final alignment.
#				(Used for -tuneup)
#	-include_all 0		In each call to select_seqs, include all
#				sequences from the HOMOLOGS set, no
#				matter how bad the scores.
#				(Used for -tuneup)

# NOT IMPLEMENTED YET:
#	-constraints 2		2 means constrain alignment on each iteration
#				to full alignment at beginning of iteration
#	-sw_select		Use Smith-Waterman instead of WU-BLAST
#				for initial selection.

# WARNING: This script is designed to use new features of SAM version 3
# It will NOT work with earlier versions, though they may
# not give error messages.
#

# Markd's notes of things to do:
#   o Script makes too much noise (mostly SAM needs to be shut up).  Problems
#     are often missed because of the noise.
#   o Should stop on all errors. This has been improved with a common proc.
#   o Need option to enable tracing of SAM commands (and arguments) that
#     are executed (currently echos all to STDERR).

# Kevin's notes of things to do:
#   o Counting nr should not be necessary---try looking it up from an index
#	(currently easiest is pcd/nrp/formatdb.log, but that is not
#	guaranteed to exist and has a poor name).
#   o Maybe eliminate use of constraints???
#   o Try to make -constraints 2 work (constraining alignment on each iteration),
#	by using new domain-map feature	Mark Diekhans added to hmmscore.

# INSTALLATION DEPENDENCIES:
#	o The first line of the file should give the path name for perl5.
#	o The perl module SamT2k.pm is loaded from the same directory containing
#         this file.
#       o It will in turn load sam-t2k.conf, which has directory
#	  information for programs and files.
#

# subroutine prototypes
sub process_command_line;
sub echo_parameters();
sub copy_initial_alignment($);
sub build_initial_alignment($);
sub check_initial_alignment();
sub get_initial_constraints($$) ;
sub select_and_train($$$$$$$$$$$);
sub select_seqs($$$$$$$$$$$$);
sub thin_a2m($$$);
sub retrain($$$);
sub realign($$$$$$$);
sub count_seqs($);
sub superfamily_params();
sub family_params();
sub tuneup_params();
sub aweight_bits($);

{   use FileHandle;
    STDOUT->autoflush(1);
    STDERR->autoflush(1);
    use English;
    use File::Basename;
    use lib dirname($PROGRAM_NAME);

    use SamT2k;

    # Set up defaults
    @search_files = ();
    $CLOSE = "";
    $HOMO = "";
    $save_all = 0;
    $no_search = 0;
    $no_sort=0;
    $force_seed = 1;
    $keep_temporary=0;
    $use_blast = 1;	# CHANGE DEFAULT TO SW WHEN WOPRKING
    $blast_prefilter  = "ncbi-blast-prefilter";
    
    $init_adpstyle = 1;
    $search_adpstyle = 1;
    $realign_adpstyle = 1;
    $final_adpstyle = 5;
    
    # adjustments to threshold for scoring against reversed sequence
    $reverse_diff = 4;

    $which_database = $NR;	# default 

    superfamily_params();
    &process_command_line;

    $use_binary = 1;	#switch to 0 if readable models needed for debugging.
    $ALLPARAM = "-fimtrans $fimtrans\\\
    	-fimstrength $fimstrength\\\
        -ins_jump_conf $trans_stiffness \\\
	-match_jump_conf $trans_stiffness \\\
	-del_jump_conf $trans_stiffness \\\
	-binary_output $use_binary\\\
	-prior_library $MIXTURE \\\
	-a2mdots 0 \\\
	-a protein \\\
        -sw 2 -jump_in_prob $jump_in -jump_out_prob $jump_out";

    $MODELFROMALIGNPARAM = "$ALLPARAM \\\
    	-aweight_method $aweight_method\\\
	-aweight_exponent $aweight_exponent";
    $BUILDPARAM= "$ALLPARAM\\\
	    -modellength 0 \\\
	    -internal_weight 0 \\\
	    -anneal_length 2 \\\
	    -many_files 1 \\\
	    -Nseq 20000\\\
	    -nsurgery $nsurgery \\\
	    -fracinsert $frac_insert \\\
	    -cutinsert $cut_insert -cutmatch $cut_match";

    $hostname = `hostname`;
    chomp $hostname;
    
    $WORKDIR  = absolute("$tmp_root_dir/target2k-$hostname-$$");
    make_directories($WORKDIR);
    chdir $WORKDIR;


    if (defined($seed_a2m))
    {    copy_initial_alignment($seed_a2m);
    }
    elsif ($CLOSE ne "")
    {    build_initial_alignment($CLOSE);
    }
    else
    {    fatal("No initial alignment.  Need to specify -seed, -homologs, or -close.\n");
    }
    check_initial_alignment();

    
    $motif_cutoff = 20.0 / $length;   # normally require 20 residues in homolog
    $motif_cutoff = 0.6 if $motif_cutoff >0.6;

    $HMMSCOREPARAM="$ALLPARAM\\\
    	    -dpstyle 0 \\\
	    -subtract_null 4 -sort 2";
	    
    if (! defined($db_size))
    {	$db_size = count_seqs($which_database);
    }
    # echoing had to wait until $db_size computed
    echo_parameters();	

    
    if ($HOMO eq "" && ! $no_search)
    {	$use_blast || die "Only know how to prefilter with blast, so far\n";
        $prefilter_string = "$prefilter_thresholds[0]";
    	for ($i=1; $i<=$#prefilter_thresholds; $i++)
	{   $prefilter_string .= ",$prefilter_thresholds[$i]";
	}
    	
	$blast_params = " -prefilter_thresholds $prefilter_string "
		. " -db $which_database -blast_max_report $blast_max_report";
	$blast_params .= " -keep_temporary" if ($keep_temporary);
	
	run_prog("$blast_prefilter -tmp_dir $WORKDIR -out prefilter -seed m0.a2m $blast_params");
	
	for ($i=0; $i<=$#prefilter_thresholds; $i++)
	{   $search_files[$i] = "prefilter_" . ($i+1) . ".fa";
	}
    }
    else
    {	for ($i=0; $i<=$#prefilter_thresholds; $i++)
   	{   $search_files[$i] = ($HOMO ne "")? $HOMO : "m0.a2m";
	}
    }
    $search_files[0] = $CLOSE if ($CLOSE ne "");

    # Build initial model using only the input sequence
    # This model is not as good as I'd like, but so far haven't found a
    # better model-building method for use with single sequences.
    # Select from the CLOSE homologs training set,
    # threshold set conservatively.

    if ($use_constraints==1) 
    {	get_initial_constraints("m0.a2m", "m0.cst");
    }
    $name = "m0";
    $num_seed_seq = count_seqs("m0.a2m");
    $possible = $search_files[0];
    
    for ($iter=1; $iter<=$num_iter; $iter++)
    {   $possible = $search_files[$iter-1];
	$num_seq = $num_seed_seq + 
		count_seqs($possible);
	# Note: it would probably be best to find one transition regularizer
	# and use it consitently
	$trans_reg= ($iter==$num_iter?  $final_trans_reg:
		    ($iter==1? $initial_trans_reg: $middle_trans_reg));

	$prev=$name;
	$name = ($save_all? $root : "m") ."_" . $iter;
	&select_and_train("$prev.a2m", "$prev.cst", "$name.a2m", "$name.cst",
                          "tmp_$iter", $possible, $num_seq, $db_size,
			  $thresholds[$iter-1], aweight_bits($iter-1),
			  $trans_reg);
    }

    if ($num_iter<1) 
    {    $num_seq = $num_seed_seq + count_seqs($search_files[0]);
    }
    select_seqs("$name.a2m", "$name.cst", 
    	"tmp_$iter", $possible,  $num_seq, $db_size,
		 $thresholds[$num_iter-1], aweight_bits($num_iter-1),
		 $final_trans_reg, $full_seq_align, $include_all, $final_adpstyle);
    
    # Have to do a copy rather than a rename for the final alignment,
    # because it may be across different file systems.
    $dest_a2m = $root . ".a2m";
    run_prog("cp tmp_$iter-a.train.seq $dest_a2m",
			"copy tmp_$iter-a.train.seq  to $dest_a2m failed");
    $dest_cst = $root . ".cst";
    if ($use_constraints) 
    {   run_prog("cp tmp_$iter-a.train.cst $dest_cst",
			"copy tmp_$iter-a.train.cst  to $dest_cst failed");
    }

    chdir "$WORKDIR/..";
    run_prog("rm -rf $WORKDIR < /dev/null",
    		"removing temporary $WORKDIR failed") if (! $keep_temporary);

    exit 0;
}

#--------------------------------------------------
# subroutines

# parse arguments

sub process_command_line
{   local($i) = 0;

    while ($i <= $#ARGV) 
    {	$_ = $ARGV[$i++];
	if (/^-seed$/)	{   $seed_a2m  = absolute( $ARGV[$i++]);
			    (-e $seed_a2m) ||
				fatal("$seed_a2m doesn't exist");
			}
	elsif (/^-out/)	{   $root = absolute($ARGV[$i++]);	}
	elsif (/^-iter/){   $num_iter	=  $ARGV[$i++];		}
	elsif (/^-jump_in/){   $jump_in	=  $ARGV[$i++];		}
	elsif (/^-jump_out/){   $jump_out	=  $ARGV[$i++];		}
	elsif (/^-fimstrength/)
			{   $fimstrength =  $ARGV[$i++];	}
	elsif (/^-fimtrans/)
			{   $fimtrans =  $ARGV[$i++];	}
	elsif (/^-all/) {   $save_all = 1;			}
	elsif (/^-init_adpstyle/)
			{   $init_adpstyle =  $ARGV[$i++];	}
	elsif (/^-search_adpstyle/)
			{   $search_adpstyle =  $ARGV[$i++];	}
	elsif (/^-realign_adpstyle/)
			{   $realign_adpstyle =  $ARGV[$i++];	}
	elsif (/^-final_adpstyle/)
			{   $final_adpstyle =  $ARGV[$i++];	}
	elsif (/^-aweight_method/)
			{   $aweight_method =  $ARGV[$i++];	}
	elsif (/^-aweight_bits/)
			{   $aweight_bitspec =  $ARGV[$i++];
			    @aweight_bitarray = split(/,/ , $aweight_bitspec);
			}
	elsif (/^-aweight_exponent/)
			{   $aweight_exponent =  $ARGV[$i++];	}
	elsif (/^-blast_max_report/)
			{   $blast_max_report =  $ARGV[$i++];	}
	elsif (/^-db_size/)
			{   $db_size = $ARGV[$i++];		}
	elsif (/^-reverse_diff/)
			{   $reverse_diff =  $ARGV[$i++];	}
	elsif (/^-thresholds/)
	{   $threshold_spec= $ARGV[$i++];
	    @thresholds = split(/,/ , $threshold_spec);
	}
	elsif (/^-prefilter_thresholds/)
	{   $prefilter_threshold_spec= $ARGV[$i++];
	    @prefilter_thresholds = split(/,/ , $prefilter_threshold_spec);
	}
	elsif (/^-tmp_dir/)
	{    $tmp_root_dir = absolute($ARGV[$i++]);
	}
	elsif (/^-homologs/)
	{    $HOMO = absolute($ARGV[$i++]);
	     (-e $HOMO) || fatal("$HOMO doesn't exist");
        }
	elsif (/^-close/)
	{    $CLOSE = absolute($ARGV[$i++]);
	     (-e $CLOSE) || fatal("$CLOSE doesn't exist");
        }
	elsif (/^-init_trans_reg/) 
	{   $initial_trans_reg = absolute_or_in_dir($ARGV[$i++], $reg_lib_dir);
	}
	elsif (/^-mid_trans_reg/) 
	{   $middle_trans_reg = absolute_or_in_dir($ARGV[$i++], $reg_lib_dir);
	}
	elsif (/^-final_trans_reg/) 
	{   $final_trans_reg = absolute_or_in_dir($ARGV[$i++], $reg_lib_dir);
	}
	elsif (/^-mixture/) 
	{   $MIXTURE = absolute_or_in_dir($ARGV[$i++], $reg_lib_dir);
	}
	
		
	elsif (/^-superfamily/)	{   superfamily_params();	}
	elsif (/^-family/)	{   family_params();		}
	elsif (/^-tuneup/)	{   tuneup_params();	}
	elsif (/^-no_search/)	{   $no_search = 1;		}
	elsif (/^-no_sort/)	{   $no_sort = 1;		}
	elsif (/^-full_seq_align/){   $full_seq_align = $ARGV[$i++];	}
	elsif (/^-include_all/)	{   $include_all = $ARGV[$i++];	}
	elsif (/^-force_seed/)	{   $force_seed = $ARGV[$i++];	}
	elsif (/^-nsurgery/)	{   $nsurgery   = $ARGV[$i++];	}
	elsif (/^-frac_insert/)	{   $frac_insert = $ARGV[$i++];	}
	elsif (/^-cut_insert/)	{   $cut_insert = $ARGV[$i++];	}
	elsif (/^-cut_match/)	{   $cut_match = $ARGV[$i++];	}
	elsif (/^-constraints/) {   $use_constraints = $ARGV[$i++];	}
	elsif (/^-ncbiblast/) {   $blast_prefilter = "ncbi-blast-prefilter";	}
	elsif (/^-wublast/) {   $blast_prefilter = "wu-blast-prefilter";	}
	elsif (/^-db/)	{   $which_database = $ARGV[$i++];
			    if ($which_database eq "nr") {$which_database=$NR;}
			    else {$which_database = absolute($which_database);}
			}
	elsif (/^-keep_temporary/)	{   $keep_temporary = $ARGV[$i++]; }
	else	{ fatal("Invalid option '$ARGV[$i-1]'");
		}
    }

    # Make certain that required parameters were provided.
    defined($root) || fatal("Must specify '-out root' with $PROGRAM_NAME");

    if ($HOMO eq "")
    {    (-e $which_database) || fatal("Can't find database file $which_database");
    }

    # verify that $num_iter>0
    if ($num_iter <=0)
    {	print STDERR "Warning: -iter specifed with number <=0, using 1\n";
	$num_iter = 1;
    }

    # Pad @thresholds if it doesn't have the right number of thresholds
    $num_thresh =  $#thresholds +1;
    if ($num_thresh  <$num_iter)
    {	print STDERR "Warning:	# thresholds ($num_thresh) < # iterations ($num_iter)\n";
    }
    for ($i=0; $i<=$#thresholds; $i++)
    {   if ($thresholds[$i] <=0 )
	{   print STDERR "$PROGRAM_NAME: Thresholds must be >0, not $thresholds[$i]\n";
	    @thresholds= @thresholds[0..$i-1];
	}
    }
    if ($#thresholds==-1)
    {   push @thresholds, -40;
    }
    while ($#thresholds < $num_iter-1)
    {   push @thresholds, ($thresholds[$#thresholds] * 0.75);
    }

    #  Make sure thare are as many prefilter_thresholds as threhsolds,
    # and that any added are at least 10 times as loose as the real thresholds.
    while ($#prefilter_thresholds < $#thresholds)
    {	push @prefilter_thresholds,
    	 (10.*$thresholds[$#prefilter_thresholds+1] > 
	 		$prefilter_thresholds[$#prefilter_thresholds])?
    		10.*$thresholds[$#prefilter_thresholds+1] 
		:$prefilter_thresholds[$#prefilter_thresholds];
    }
    
    
    $CLOSE = $HOMO if ($CLOSE eq "");
    defined($seed_a2m) 
	|| defined($CLOSE) 
	|| defined($HOMO) 
    	|| fatal("Must specify '-seed file.a2m', '-close file.seqs', or '-homologs file.seqs' with $PROGRAM_NAME");

    if ($force_seed>0 && $nsurgery >0)
    {    print STDERR "-force_seed $force_seed is incompatible with "
    		. "-nsurgery $nsurgery\n"
		. "Setting -nsurgery 0\n";
    	$nsurgery = 0;
    }

}


sub echo_parameters()
{    # Echo the parameters to stderr
    print STDERR "$PROGRAM_NAME: E-value thresholds used are\n";
    for ($i=0; $i<$num_iter; $i++)
    {	print STDERR "prefilter: $prefilter_thresholds[$i],\t"
    			."actual: $thresholds[$i]\n";
    }

    print STDERR "$PROGRAM_NAME: aweight_bits are ", aweight_bits(0);
    for ($i=1; $i<$num_iter; $i++)
    {	print STDERR ", " , aweight_bits($i);
    }
    print STDERR "\n";

    print STDERR "$PROGRAM_NAME: db_size= $db_size reverse_diff= $reverse_diff\n";
    print STDERR "$PROGRAM_NAME: blast_prefilter: $blast_prefilter\n";
    if ($HOMO eq "")
    {	print STDERR "$PROGRAM_NAME: blast_max_report= $blast_max_report\n";
    	print STDERR "$PROGRAM_NAME: db=$which_database\n";
    }
    else
    {	print STDERR "$PROGRAM_NAME: homologs=$HOMO close=$CLOSE\n";
    }
    print STDERR "$PROGRAM_NAME: aweight_method=$aweight_method aweight_exponent=$aweight_exponent\n";
    print STDERR "$PROGRAM_NAME: force_seed=$force_seed constraints=$use_constraints\n";
    print STDERR "$PROGRAM_NAME: jump_in=$jump_in jump_out=$jump_out fimstrength=$fimstrength fimtrans=$fimtrans\n";
    if ($nsurgery>0)
    {	print STDERR "$PROGRAM_NAME: nsurgery=$nsurgery cutinsert=$cut_insert cutmatch=$cut_match fracinsert=$frac_insert\n";
    }
    return;
}

# Look up alignment and copy to init.a2m
sub copy_initial_alignment($)
{
    my($from_seed_a2m) = @_;
    $BASE_A2M = ($from_seed_a2m =~ /(\S+)[.]gz$/ )? $1 : $from_seed_a2m;
    if (-e $BASE_A2M)
    {    run_prog("cp $BASE_A2M m0.a2m", 
    		"Copying initial alignment to m0.a2m failed");
    }
    elsif (-e "$BASE_A2M.gz")
    {   # Use gunzip to avoid picking up the wrong zcat
        run_prog("gunzip -c $BASE_A2M.gz > m0.a2m",
    		"uncompress to m0.a2m failed");
    }
    else
    {   fatal("Error: can't find seed alignment \"$BASE_A2M\" or \"$BASE_A2M.gz\"");
    }
}


# Build initial alignment from a set of unaligned homologs
# leaving the result in m0.a2m
sub build_initial_alignment($)
{
    my ($unaligned_file) = @_;
    # build a model by attempting to align the unaligned sequences
    run_prog("buildmodel m0 -train $unaligned_file \\\
    	$BUILDPARAM \\\
		-initial_noise 0.05 -anneal_noise 0.01 \\\
		-Nmodels 7 -sequence_models 0.5\\\
		-insert $initial_trans_reg",
	"buildmodel of unaligned sequences failed");
    
    # align all the sequences, no matter how they score.
    run_prog("hmmscore m0 -i m0.mod -db $unaligned_file \\\
	$ALLPARAM \\\
	-adpstyle $init_adpstyle \\\
	-select_score 8 -select_align 8",
		"hmmscore failed");


}


# check the initial alignment and get the number of alignment columns
# (returned as $length)
sub check_initial_alignment()
{    
    run_prog("checkseq foo -alphabet protein -db m0.a2m > init.check",
             "checkseq failed on initial alignment");

    # TO DO:  if m0.a2m fails checks, then build initial alignment
    #	from sequences of init.a2m
    
    open(CHECK, "<init.check");
    while(<CHECK>)
    {	if (/^Alignment:\s+no/)
	{    fatal("m0.a2m not recognized as an alignment");
	}
	if (/^AlignColumns:\s+(\d+)/)
	{   $length = $1;
	    if ($length <= 1)
	    {   fatal("m0.a2m has $length alignment columns\n" .
                      "A2M format requires upper-case for aligned residues");
	    }
	}
    };
    close CHECK;
}


#
# Create initial constraints for initial alignment.  Although this
# could be done by the first select_and_train, doing it upfront
# keeps that code cleaner.
#   o inita2m - Input a2m
#   o initcst - Output constraints.
#
sub get_initial_constraints($$) 
{   my ($align_in, $constraints_out) = @_;
    run_prog("modelfromalign unused -alignfile $align_in \\\
             -constraints_from_align 1 -constraints_out $constraints_out",
	     "modelfromalign for initial constraints failed");
}


# Do one round of selection and training, starting and ending with an alignment
#	o name of initial alignment
#       o name of initial alignment's constraints file
#	o name of final alignment
#       o name of final alignment's constraints file
#	o base name for temporaries
#	o set of possible homologs to search
#	o max number of sequences in possible+initial
#	o number of sequences in database to assume for E-value computations
#	o threshold to use
#	o weighting bits to save
#	o transition regularizer to use
sub select_and_train($$$$$$$$$$$)
{    my ($align_in, $constraints_in, $align_out, $constraints_out,
	$base, $homologs, $num_homolog_seqs,
	$num_db_seqs, $e_thresh, $aweight_bits, $trans_reg) = @_;

    select_seqs($align_in, $constraints_in,
	$base, $homologs, $num_homolog_seqs, $num_db_seqs, $e_thresh,
		$aweight_bits, $trans_reg, $full_seq_align, $include_all, $search_adpstyle);
    # bug fixed 20 march 2002:
    #	Was previously always using mult-domain, ignoring $full_seq_align.
    
    thin_a2m("$base-a.train.seq", $nll_thresh, "$base-a.train.seq");
    retrain($base, $aweight_bits, $trans_reg);   
    realign($align_out, $constraints_out, $base, $db_size, 
    	$e_thresh, $trans_reg, $realign_adpstyle);   
}

# build_hmmscore_cmd is intended to be called only from within
# select_sequences---parameters are identical
sub build_hmmscore_cmd($$$$$$$$$$$$)
{    my ($align_in, $constraints_in, 
	$base, $homologs, $num_homolog_seqs, $num_db_seqs,
	$e_thresh, $aweight_bits, $trans_reg, $full_seq_only,
	$include_all, $adpstyle) = @_;
    
#    print STDERR "DEBUG: include_all=$include_all\n";
    
    my $selection_param = $include_all? 8 : 4;

    $cmd = "hmmscore $base-a -i $base-a.mod  $HMMSCOREPARAM \\\
	-db_size $num_db_seqs \\\
    	-adpstyle $adpstyle \\\
	-align_short 0	\\\
	-simple_threshold $simple_thresh -Emax $e_thresh \\\
	-mdEmax $e_thresh\\\
        -select_score $selection_param \\\
	-select_seq $selection_param ";

    $cmd .=  $full_seq_only? " -select_align $selection_param"
    		: " -select_mdalign $selection_param -select_md $selection_param -Motifcutoff $motif_cutoff";
			# ??? tickle bugs if $selection_param

	
    if (! $force_seed && defined($seed_a2m))
    {	$cmd .= " -db $seed_a2m";
    }

    if ($homologs ne "")
    {	$cmd .= " -db $homologs";
    }
    if ($use_constraints) 
    {	$cmd .= " -constraints $base-a.cst";
	$cmd .= " -constraints_out $base-a-md.cst";
    }
    return $cmd;
}

# select domains from $homologs (and maybe $align_in)
# Leave alignment in file $base-a.trans.seq
# 	and constraints in $base-a.trans.cst
sub select_seqs($$$$$$$$$$$$)
{    my ($align_in, $constraints_in, 
	$base, $homologs, $num_homolog_seqs, $num_db_seqs,
	$e_thresh, $aweight_bits, $trans_reg, $full_seq_only,
	$include_all, $adpstyle) = @_;

    $db_over_expected = ($num_db_seqs +0.0)/ $e_thresh ;
    # compute thresholds from E values.
    # NOTE: thresholds are global, since they are used in select_and_train.
    $nll_thresh = $db_over_expected<=1.0000001? 16.: -log($db_over_expected-1);
    $simple_thresh = $nll_thresh - $reverse_diff;

    $align_for_model = $align_in;
#    thin_a2m($align_in, $nll_thresh, "$base-a-thin_init.a2m");
#    $align_for_model = "$base-a-thin_init.a2m";
    
    # build inital model
    $cmd = "modelfromalign $base-a -alignfile $align_for_model \\\
		-insert $trans_reg \\\
		-aweight_bits $aweight_bits\\\
		$MODELFROMALIGNPARAM";
    if ($use_constraints !=0)
    {   if (defined($constraints_in))
	{   $cmd .= " -constraints $constraints_in";
	}
	else
	{   $cmd .= " -constraints_from_align 1";
	}
	$cmd .= " -constraints_out $base-a.cst";
    }
    run_prog($cmd, "modelfromalign for initial model failed");


    if (($force_seed || !defined($seed_a2m)) && $homologs ne "")
    {   if (-z $homologs)
	{   # oops no known homologs in $homologs, nothing to look for.
	    run_prog("cp $homologs  $base-a-sort.a2m",
		"empty copy of $homologs failed");
	    run_prog("cp $seed_a2m $base-a.train.seq",
	    	"copying seed to next iteration failed");
	    run_prog("cp $seed_a2m $base-a.sel",
	    	"copying seed to sequence set for realignment failed");
	    if ($use_constraints)
	    {   run_prog("cp m0.cst $base-a.train.cst",
			"copying seed constraints failed");
            }
	    return;
	}
    }

    # select from possible homologs (using simple null and reverse sequence)
    # then reselect only the relevant parts
    #	(but add some extra residues for context)
    # Keep both the full sequences and the selected parts, but not
    # alignment of full sequences, unless $full_seq_only is set.
    $cmd=build_hmmscore_cmd($align_in, $constraints_in, 
	$base, $homologs, $num_homolog_seqs, $num_db_seqs,
	$e_thresh, $aweight_bits, $trans_reg, $full_seq_only,
	$include_all, $adpstyle);
    run_prog($cmd, "hmmscore to select sequences failed");

    $selected_parts = $full_seq_only?  "$base-a.a2m" : "$base-a.mult";
    $constraints =  $full_seq_only?  "$base-a.cst" : "$base-a-md.cst";
    $costs = $full_seq_only?  "$base-a.dist" : "$base-a.mstat";
    
    # $base-a.sel has the selected homologs (unaligned)
    # $selected_parts has parts of the selected homologs to train on
    #		(aligned)
    # $constraints has the inherited constraints on $selected_parts
    #

    if (-z "$selected_parts" && !$force_seed)
    {    # rats, got an empty search result and not forcing the seed.
	 # Select everything in this emergency!
	 $cmd=build_hmmscore_cmd($align_in, $constraints_in, 
		$base, $homologs, $num_homolog_seqs, $num_db_seqs,
		$e_thresh, $aweight_bits, $trans_reg, $full_seq_only,
		1, $adpstyle);
    	 run_prog($cmd, "hmmscore to select all sequences failed");
    }
    
    if (-z "$selected_parts" || $no_sort)
    {   # don't call sortseq
    	run_prog("cp $selected_parts  $base-a-sort.a2m",
		"copying $selected_parts to $base-a-sort.a2m failed");
    }
    else
    {   run_prog("sortseq $base-a-sort -alignfile $selected_parts \\\
    		-NLLfile $costs -a2mdots 0", 
		"sortseq failed");
    }
    # $base-a-sort.a2m has parts of the selected homologs to train on
    #		sorted best first (unless $no_sort set)

    # add the initial alignment to the training set, to make sure it doesn't
    # get lost along the way
    if ($force_seed)
    {	run_prog("cat m0.a2m $base-a-sort.a2m > $base-a.train.seq",
    		"adding seed to $base-a-sort.a2m failed");
        if ($use_constraints)
	{   run_prog("cat m0.cst $constraints > $base-a.train.cst",
		"adding seed constraints failed");
        }
    }
    else
    {	rename("$base-a-sort.a2m", "$base-a.train.seq");
        if ($use_constraints)
	{   rename("$constraints", "$base-a.train.cst");
        }
    }
}

# thin the specified alignment, removing duplicate (and near_duplicate) ids
sub thin_a2m($$$)
{   my ($thick_file, $nll_thresh, $thin_file) = @_;
    # Set fraction identical heuristically based on $nll_thresh and
    #	alignment size
    # As $nll_thresh->-inf, $frac_id->1.0,  as $nll_thresh->inf, $frac_id->$min_frac_id
    # As $num_seqs -> inf, fraction identical->$min_frac_id
    $num_seqs = count_seqs($thick_file);
    ( $num_seqs > 0 ) 
    	|| fatal("No sequences in $thick_file in thin_a2m\n" );
    $min_frac_id = 0.2;
    $frac_id =  $min_frac_id + (1.-$min_frac_id)*
    	(1.- 1./(1+exp(-10.0*$nll_thresh/$num_seqs)));
    print STDERR "# num_seqs=$num_seqs nll_thresh=$nll_thresh frac_id=$frac_id\n";

    run_prog("uniqueseq thinned_tmp -alignfile $thick_file -a protein \\\
        -percent_id $frac_id -a2mdots 0","uniqueseq failed");
    rename "thinned_tmp.a2m", $thin_file;
}



# Build a new model $base-b.mod, starting from $base-a.mod
# using training sequences in $base-1.train.seq.
sub retrain($$$)
{    my ($base, $aweight_bits, $trans_reg) = @_;
    
    if ($weighted_build)
    {	# get weights from initial alignment of training sequences
    	run_prog("modelfromalign junk_for_weights \\\
    		-alignfile $base-a.train.seq \\\
		-aweight_bits $aweight_bits\\\
    		$MODELFROMALIGNPARAM",
		"modelfromalign to get sequence weights failed");
	rename("junk_for_weights.weightoutput", "$base-a.w");
    }

    $cmd = "buildmodel $base-b -insert $base-a.mod\\\
	 $BUILDPARAM -initial_noise 0.0001 -anneal_noise 0.000001 \\\
		-train $base-a.train.seq \\\
		-nmodels 1\\\
		-insert $trans_reg";
    if 	($weighted_build)
    { 	$cmd .= " -sequence_weights $base-a.w";
    }
    if ($use_constraints)
    {	$cmd .= " -constraints $base-a.train.cst";
    }
    run_prog($cmd, "buildmodel for retraining failed");
}

# Build a new alignment in $align_out, using $base-b.mod,
# using sequences from $base-a.sel
sub realign($$$$$$$)
{    my ($align_out, $constraints_out, $base, $num_db_seqs,
	$e_thresh, $trans_reg, $adpstyle) = @_;
    
    #
    # Reselect interesting pieces out of the selection set.
    # Note: there is a slight inefficiency here, since $base-a.sel
    # may overlap signficiantly with m0.a2m.  Unfortunately,
    # there may be some important sequences in one and not the other,
    # and I haven't figured out exactly when I can know that they are
    # redundant.
    #
    $dbs = $include_all? "-db $search_files[0]" 
    		: "-db m0.a2m -db $base-a.sel";
    $selection_param = $include_all? 8 : 4;

    
    $cmd = "hmmscore $base-b -i $base-b.mod $HMMSCOREPARAM \\\
	-db_size $num_db_seqs \\\
	$dbs \\\
	-select_score $selection_param \\\
	-select_seq $selection_param \\\
	-simple_threshold $simple_thresh -Emax $e_thresh \\\
	-select_mdalign 4 -select_md 4 \\\
	-adpstyle $adpstyle \\\
	-Motifcutoff $motif_cutoff \\\
	-align_short 5	-mdEmax $e_thresh";
	# tickle bugs if $selection_param for -select_md    
    if ($use_constraints)
    {	$cmd .= " -constraints $base-a.train.cst";
        $cmd .= " -constraints_out $constraints_out";
    }

    run_prog($cmd,"hmmscore to realign sequences failed");

    #
    # Return output files.
    #
    $final_align = "$base-b.mult";
    if ("$final_align" ne $align_out)
    {	run_prog("cp -f $final_align $align_out", 
    		"copying realigned sequences failed");
    }

    if (-z $align_out)
    {	print STDERR "*** Warning: Alignment $align_out empty, copying initial alignment\n";
	run_prog("cp -f m0.a2m $align_out", "copying seed failed");
        if ($use_constraints) 
	{   run_prog("cp -f m0.cst $constraints_out", "copying seed constraints failed");
        }
    }
}

#
#  count the number of sequences in a FASTA file
#
sub count_seqs($)
{   my ($file_name) = @_;
    return 0 if (!defined($file_name) || $file_name eq "");
    print STDERR "# About to count sequences in $file_name\n";
    $wc_ret = `grep '^>' $file_name | wc`;
    $wc_ret =~ /\s*(\d+)/;
    return $1;
}


# Set parameters appopriately for a superfamily alignment starting
# from a single sequence
sub superfamily_params()
{
    $num_iter=4;
    # thresholds for E-value for full NR size
    @thresholds= (0.0001, 0.0002, 0.001, 0.005);	# for HMMSCORE
    @prefilter_thresholds = (0.01, 1.0, 10.0, 400.0);	
    
    $blast_max_report= 20000;
    
    $jump_in = 0.2;
    $jump_out = 1.0;
    $fimstrength= 1.0;
    $fimtrans = -1.0;
    
    $use_constraints = 1;
    $force_seed=1;
    
    $weighted_build = 1;
    $aweight_method = 1;
    @aweight_bitarray = (0.8,0.7,0.6,0.5);
    $aweight_exponent = 10.0;
    $aweight_var_bits=0;
    
    $trans_stiffness = 1.0;
    
    $nsurgery=0;
    $frac_insert = 0.2;
    $cut_insert = 0.25;
    $cut_match = 0.2;
    
    $MIXTURE = "$reg_lib_dir/recode3.20comp";
    $initial_trans_reg = "$reg_lib_dir/gap1.5.regularizer";
    $middle_trans_reg = "$reg_lib_dir/stiff-gap5.regularizer";
    $final_trans_reg = "$reg_lib_dir/fssp-trained.regularizer";
    
    $full_seq_align=0;
    $include_all =0;
}

# Set parameters appopriately for a family alignment starting
# from a single sequence.
# WARNING: THIS IS ONLY A FIRST GUESS
sub family_params()
{
    $num_iter=3;
    # thresholds for E-value
    @thresholds= (1.e-8, 1.e-6, 1.e-5);
    @prefilter_thresholds=(1.e-10, 1.e-3, 10.);
    
    $blast_max_report= 10000;

    $jump_in = 0.2;
    $jump_out = 1.0;
    $fimstrength= 1.0;
    $fimtrans = -1.0;

    $use_constraints = 1;
    $force_seed=1;
    
    $weighted_build = 1;
    $aweight_method = 1;
    @aweight_bitarray = (0.8,0.7,0.5);
    $aweight_exponent = 10.;
    $aweight_var_bits=0;
    
    $trans_stiffness = 1.0;

    $nsurgery=0;
    $frac_insert = 0.2;
    $cut_insert = 0.25;
    $cut_match = 0.2;
    
    $MIXTURE = "$reg_lib_dir/recode3.20comp";
    $initial_trans_reg = "$reg_lib_dir/gap1.5.regularizer";
    $middle_trans_reg = "$reg_lib_dir/stiff-gap5.regularizer";
    $final_trans_reg = "$reg_lib_dir/fssp-trained.regularizer";
    
    $full_seq_align=0;
    $include_all =0;
}


# Set parameters appopriately for tuning up an initial alignment.
sub tuneup_params()
{
    $num_iter=3;
    # thresholds for E-value
    @thresholds= (0.0001, 0.001, 0.01);	# may be irrelevant, with $include_all
    	# except for $nll_thresh used in thinning
    $db_size=10;	# should probably be size of seed+homologs.
    
    $no_search=1;	$no_sort=1;
    $jump_in = 0.2;
    $jump_out = 1.0;

    $fimstrength= 1.0;
    $fimtrans = -1.0;
    $use_constraints = 0;
    $force_seed=0;
    
    $weighted_build = 1;
    $aweight_method = 1;
    $aweight_exponent = 10;
    $aweight_var_bits=0;
    @aweight_bitarray = (0.5,0.5,0.5);	
    
    $trans_stiffness = 1.0;
    
    $nsurgery=10;
    $frac_insert = 0.3;
    $cut_insert = 0.3;
    $cut_match = 0.2;
    
    $MIXTURE = "$reg_lib_dir/recode3.20comp";
    $initial_trans_reg = "$reg_lib_dir/gap1.5.regularizer";
    $middle_trans_reg = "$reg_lib_dir/stiff-gap5.regularizer";
    $final_trans_reg = "$reg_lib_dir/fssp-trained.regularizer";
    
    $full_seq_align=1;
    $include_all =1;
}


# Return target bits/column for iteration 0..n-1
sub aweight_bits($)
{    my($iter) = @_;
    my ($aweight_bits, $aweight_var_bits);
    $aweight_bits =  defined($aweight_bitarray[$iter])?
    	$aweight_bitarray[$iter] :
    	$aweight_bitarray[$#aweight_bitarray];	
    if ($aweight_bits =~ /^[vV][aA][rR](\S+)/)
    {   $aweight_var_bits = $1;
        $pred_len_align = 140.0* (1.0-exp(-0.008*$length));
	$pred_len_align = $length if $pred_len_align > $length;
	# Attempting a bug fix:
	# The number of bits to save got too large for very
	#    short alignments, so bound length from below.
	$pred_len_align = 50.0 if $pred_len_align < 50.0;
	$aweight_bits = (0.0+$aweight_var_bits) / $pred_len_align;
    }
    return $aweight_bits;
    
}
