* Fix perl interpreter paths.
* Support parallel build of binaries.
* Enhance verbosity in csh scripts for bug hunting.
* Add blast processors argument SIFT_for_submitting_fasta_seq.csh .
* Trailing whitespace removed from lines.
Index: sift/bin/perlscripts/get_BLINK_seq.pl
===================================================================
--- sift.orig/bin/perlscripts/get_BLINK_seq.pl
+++ sift/bin/perlscripts/get_BLINK_seq.pl
@@ -1,4 +1,4 @@
-#!/usr/local/bin/perl 
+#!/usr/bin/perl
 #use strict;
 use LWP::Simple;
 use LWP::UserAgent;
Index: sift/bin/seqs_chosen_via_median_info.csh
===================================================================
--- sift.orig/bin/seqs_chosen_via_median_info.csh
+++ sift/bin/seqs_chosen_via_median_info.csh
@@ -1,7 +1,8 @@
-#!/bin/csh
+#!/bin/csh -e
 # seqs_chosen_via_median_info.csh
 # Arg1 = protein sequence in fasta format
 # Arg2 = blastpgp database
+# Arg3 = blastpgp processors (-a)
 
 # This program is licensed to you under the Fred
 # Hutchinos Cancer Research Center (FHCRC)
Index: sift/bin/SIFT_for_submitting_NCBI_gi_id.csh
===================================================================
--- sift.orig/bin/SIFT_for_submitting_NCBI_gi_id.csh
+++ sift/bin/SIFT_for_submitting_NCBI_gi_id.csh
@@ -1,4 +1,4 @@
-#!/bin/csh
+#!/bin/csh -e
 #		SIFT.csh
 # This program is licensed to you under the Fred
 # Hutchinos Cancer Research Center (FHCRC)
Index: sift/src/Alignment.c
===================================================================
--- sift.orig/src/Alignment.c
+++ sift/src/Alignment.c
@@ -17,7 +17,7 @@ query sequence is moved to be the first
 #define _ALIGNMENT_C_
 #include "PN_convert.c"
 
-#define MAXSEQ 400
+#define MAXSEQ 20000
 
 #include "Alignment.h"
 
@@ -478,6 +478,7 @@ int flat_master_slave (FILE* fp, Sequenc
 		fprintf(errorfp, "ERROR! Why haven't %s and %s been freed?\n", prev_name,
 		prev_residues);
 		fprintf (errorfp, "Current %s %s\n", name, residues);
+		fprintf (stderr, "Current %s %s\n", name, residues);
 		exit (-1);
 	 }
 	} /* end of if (line[0] != ' ') */ 
@@ -752,7 +753,8 @@ read_psiblast_header_until_first (FILE*
         }
 	if (feof (fp)) {
 		fclose (fp);
-			fprintf (errorfp, "ERROR! PSI-BLAST found no hits. Program terminating.\n");
+		fprintf (errorfp, "ERROR! PSI-BLAST found no hits. Program terminating.\n");
+		fprintf (stderr, "%s:%d ERROR! PSI-BLAST found no hits. Program terminating.\n", __FILE__, __LINE__ );
                 exit (-1);
 	}
 
@@ -792,7 +794,8 @@ read_psiblast_header_until_first_no_erro
                 fclose (fp);
                 if (return_error) {
                         fprintf (errorfp, "ERROR! PSI-BLAST found no hits. Program terminating.\n");
-                exit (-1);
+                        fprintf (stderr, "%s:%d ERROR! PSI-BLAST found no hits. Program terminating.\n", __FILE__, __LINE__ );
+                	exit (-1);
                 } else {
                         return (-1);
                 }
@@ -827,6 +830,7 @@ printf ("finished reading psiblast file\
 	if (feof (fp)) {
 		fclose (fp);
 		fprintf (errorfp, "ERROR! PSI-BLAST found no hits. Program terminating.\n");
+		fprintf (stderr, "%s:%d ERROR! PSI-BLAST found no hits. Program terminating.\n", __FILE__, __LINE__ );
 		exit (-1);
 	}
 	sscanf (line, "Results from round %d", &iteration);
@@ -1213,6 +1217,7 @@ get_index_from_seqs (char name[], Sequen
 		}
 	}
 	fprintf (errorfp, "EXITING -- couldn't find %s in get_index_from_seqs\n", name);
+	fprintf (stderr, "EXITING -- couldn't find %s in get_index_from_seqs\n", name);
 	exit (-1);
 }
 
@@ -1352,8 +1357,9 @@ psiblast_pairwise (FILE* fp, Sequence *
 		sscanf (strptr, "Expect = %lf", &evalue);
 		} else { 
 			fprintf (errorfp, "Unable to read e-value, parsing incorrect");
-                exit (-1);
-        }
+			fprintf (stderr, "Unable to read e-value, parsing incorrect");
+                	exit (-1);
+        	}
 /* printf ("evalue is %lf\n", evalue); */ 
 /* end get evalue */
 			
@@ -1439,6 +1445,7 @@ nonoverlapping (Sequence* seq1, Sequence
 	assert (seq1->length == seq2->length);
 	if (strcmp (seq1->name, seq2->name ) != 0) {
 		fprintf (errorfp, "seeing whether %s and %s overlap, should be same name\n", seq1->name, seq2->name);
+		fprintf (stderr, "seeing whether %s and %s overlap, should be same name\n", seq1->name, seq2->name);
 		exit (-1);
 	}
 	for (i = 0; i < seq1->length; i++) {
@@ -1543,7 +1550,7 @@ read_4lines_of_pairwise (FILE* fp, Seque
 	}
 
 	if (strstr (line, "Query") == NULL) {
-		printf ("%s line should have Query:\n", line);
+		fprintf (stderr, "%s line should have Query:\n", line);
 		exit (-1);
 	}
         strptr = strtok (line, " \t\r\n"); /* this should be Query */
@@ -2052,10 +2059,10 @@ printf ("entered if statment\n");
    }
 
 
- printf ("Clumped %d sequences into %d clusters\n", nseq, i2);
+ printf("Clumped %d sequences into %d clusters\n", nseq, i2);
  /* keep this print because get # of clustersin awk statement */
   free(pairs);
-printf ("exiting clump\n");
+fprintf (stderr, "%s:%d exiting clump\n", __FILE__, __LINE__);
  return no_of_clumps;
 }  /* end of cluster */
 
Index: sift/src/blocklist.c
===================================================================
--- sift.orig/src/blocklist.c
+++ sift/src/blocklist.c
@@ -9,7 +9,7 @@
 #include <string.h>
 #include "blocklist.h"
 
-#define MAXSEQ 400
+#define MAXSEQ 20000
 
 /*======================================================================
 routines for a list of blocks
Index: sift/src/choose_seqs_via_psiblastseedmedian.c
===================================================================
--- sift.orig/src/choose_seqs_via_psiblastseedmedian.c
+++ sift/src/choose_seqs_via_psiblastseedmedian.c
@@ -1,6 +1,6 @@
 /* (C) Copyright 2000, Fred Hutchinson Cancer Research Center */
 /* Use, modification or distribution of these programs is subject to */
-/* the terms of the non-commercial licensing agreement in license.h. */           
+/* the terms of the non-commercial licensing agreement in license.h. */
 /* choose_seqs_via_psiblast.c */
 /* 02-14-01 use median threshold
 05-18-01, changed to add PSIBLAST_TURN sequences at once, if threshold has not
@@ -8,7 +8,7 @@
 	 recalculate median threshold for the 1 .. PSIBLAST_TURN sequences
 	 that have been added.
 	Should reduce time for calculation of median threshold O (length * seq)
-	  by PSIBLAST_TURN x fold	 
+	  by PSIBLAST_TURN x fold
 */
 #define EXTERN
 #include "blocksprogs.h"
@@ -20,7 +20,7 @@
 #include "Matrix_Info.c"
 #include "Protdist.c"
 #include "stringhash.c"
-#include "Psiblast.c" 
+#include "Psiblast.c"
 #include "PN_blocks.c"
 #include "Clumping.c"
 
@@ -29,9 +29,9 @@
 
 void getargs (int argc, char* argv[],
 		 char query_seq_file[LARGE_BUFF_LENGTH],
-		FILE** seqfp, 
+		FILE** seqfp,
 		FILE** outfp,
-		double* clumping_threshold, 
+		double* clumping_threshold,
 		char pid[SMALL_BUFF_LENGTH],
 		double* median_threshold);
 
@@ -53,10 +53,10 @@ main (int argc, char* argv[])
 	FILE* outfp; FILE* temp_chk_fp;
 	FILE* seqfp;
 	int nseqs, db_type, seq_type;
-	int i, k, index, nseqs_written; double dtemp;
+	int i, k, index = 0, nseqs_written; double dtemp;
 	double clumping_threshold;
 	double info_median;
-	Block* workingblock; 
+	Block* workingblock;
 	HashTable seqnamehash;
 	char seqname[KEY_WIDTH];
 	char pid[SMALL_BUFF_LENGTH];
@@ -226,15 +226,15 @@ error file */
 /*		print_block_ids (workingblock, outfp, TRUE); */
 		/* should free memeory, clean up later */
 /*		fclose (outfp); fclose (errorfp); fclose (logfp);
-		exit (0);*/ 
+		exit (0);*/
 	/* exit 0 else will generate an error in shell that calls this program*/
 	}
 	/* now of the last PSIBLAST_TURN sequences added, find out which one gives
-	   approaches the right median */	
+	   approaches the right median */
 	num_seqs_removed = 0;
 	while (median < median_threshold) {
 		workingblock->num_sequences--;
-		median = calculate_median_information_of_block (workingblock, 
+		median = calculate_median_information_of_block (workingblock,
 							FALSE, FALSE);
 		num_seqs_removed++;
 	}
@@ -243,30 +243,30 @@ error file */
 
 	fprintf (logfp, "********FINAL BLOCK ****** ");
 	fprintf (logfp, "last block : median %.3f with %d sequences\n", median,
-			workingblock->num_sequences); 
-/*	output_block_s (workingblock, outfp, FLOAT_OUTPUT); */ 
+			workingblock->num_sequences);
+/*	output_block_s (workingblock, outfp, FLOAT_OUTPUT); */
 
 
-/*	print_block_sequences (workingblock, outfp); */ 
-	print_block_ids (workingblock, outfp, TRUE); 
+/*	print_block_sequences (workingblock, outfp); */
+	print_block_ids (workingblock, outfp, TRUE);
 
-	workingblock->num_sequences+= num_seqs_removed; 
+	workingblock->num_sequences+= num_seqs_removed;
 				/* increase it back to deallocate
-					memory properly */	
+					memory properly */
 	free_seqs (seqs, nseqs);
 	free_block (workingblock);
 	DestroyTable (seqnamehash);
 	fprintf (logfp, "SUCCESSFUL\n");
-	fclose (outfp);	
+	fclose (outfp);
 	fclose (errorfp);
-	fclose (logfp);	
+	fclose (logfp);
 	rm_file (errorfilename);
 	exit (0);
 } /* end of main */
 
 /* return the number of sequences written */
 int
-write_seqs_not_in_hash (char query_database_file[LARGE_BUFF_LENGTH], 
+write_seqs_not_in_hash (char query_database_file[LARGE_BUFF_LENGTH],
 			 Sequence* seqs[MAXSEQ], int nseqs,
                                          HashTable hash)
 {
@@ -277,7 +277,7 @@ write_seqs_not_in_hash (char query_datab
 
 	nseqs_written = 0;
 	if ( (outfp = fopen (query_database_file, "w")) == NULL) {
-		fprintf (errorfp, "couldn't write to %s \n", 
+		fprintf (errorfp, "couldn't write to %s \n",
 			query_database_file);
 		exit (-1);
 	}
@@ -291,15 +291,15 @@ write_seqs_not_in_hash (char query_datab
 
 	fclose (outfp);
 	return nseqs_written;
- 
+
 } /* end write_seqs_not_in_block */
-		
+
 void
-getargs (int argc, char* argv[], 
+getargs (int argc, char* argv[],
 	char query_seq_file[LARGE_BUFF_LENGTH],
-	FILE** seqfp, 
-	FILE** outfp, 
-	double* clumping_threshold, 
+	FILE** seqfp,
+	FILE** outfp,
+	double* clumping_threshold,
 	char pid[SMALL_BUFF_LENGTH],
 	double* median_threshold)
 {
@@ -359,7 +359,7 @@ getargs (int argc, char* argv[],
 	}
 
 	if (argc > 4) *clumping_threshold = atof (argv[4]);
-	else 
+	else
 	{
 		printf ("Enter clumping threshold\n");
 		scanf ("%lf\n", clumping_threshold);
@@ -381,10 +381,10 @@ getargs (int argc, char* argv[],
 	}
 
 
-} /* end of getargs */ 
+} /* end of getargs */
 
 void
-write_seq_to_file (char filename[LARGE_BUFF_LENGTH], Sequence* seq) 
+write_seq_to_file (char filename[LARGE_BUFF_LENGTH], Sequence* seq)
 {
 	FILE* fp;
 
Index: sift/src/clump_output_alignedseq.c
===================================================================
--- sift.orig/src/clump_output_alignedseq.c
+++ sift/src/clump_output_alignedseq.c
@@ -25,7 +25,7 @@ option 1 -- print out sequences that hav
 #include "Clumping.c"
 
 #define MINSEQ 2 
-#define MAXSEQ 400
+#define MAXSEQ 20000
 #define INDEX(n, col, row) (col*n - (col*(col+3))/2 - 1 + row)
 #define MIN_CLUSTERS 3
 
Index: sift/src/info_on_seqs.c
===================================================================
--- sift.orig/src/info_on_seqs.c
+++ sift/src/info_on_seqs.c
@@ -11,7 +11,7 @@
 #include "Psiblast.c"
 /* #include "Information.c */
 
-#define MAXSEQ 400 /* Maximum number of sequences */
+#define MAXSEQ 20000 /* Maximum number of sequences */
 #define LINE_LEN 800
 #define FALSE 0
 #define TRUE 1
@@ -200,10 +200,10 @@ output_predictions (AAnode* polymorph_da
 	int substitution, original_aa;
 	int substitution_exists;
 	int pos;
-	double median;
+	double median = 0;
 	AAnode current;
-	Block* block_with_seqs_at_pos;
-	double* info_array;
+	Block* block_with_seqs_at_pos = NULL;
+	double* info_array = NULL;
 	int block_constructed;
 	
 	substitution_exists = FALSE;
@@ -339,7 +339,7 @@ comments_on_info(Block* oldblock, int po
        fprintf (outfp, "%d\n", residues_stored[pos]);
 */
 	/* the array will be sorted in call for median_of_array*/
-printf ("looking for median\n");
+fprintf( stderr, "looking for median\n");
 	median = median_of_array (info_array, prot_length);
        printf ("found mediant %.2f\n", median);
 	 fprintf (outfp, "%.3f\t%d\t%d\n",
@@ -354,7 +354,7 @@ printf ("looking for median\n");
 	assert (oldblock->sequences[0].sequence != NULL);
 	assert (oldblock->sequences != NULL);
 	assert (oldblock != NULL);
-printf ("finished comments oninfo\n");
+fprintf( stderr, "finished comments oninfo\n");
 }
 
 double
Index: sift/src/Makefile
===================================================================
--- /dev/null
+++ sift/src/Makefile
@@ -0,0 +1,45 @@
+include-prefix := $(prefix)
+lib-prefix := $(prefix)
+
+BIN :=	choose_seqs_via_psiblastseedmedian \
+	clump_output_alignedseq \
+	consensus_to_seq \
+	info_on_seqs \
+	psiblast_res_to_fasta_dbpairwise \
+	seqs_from_psiblast_res
+
+CC := gcc
+CFLAGS := -I$(include-prefix)/include/blimps -O2 -D__MAKE_PREFIX__=\"$(prefix)\" $(CFLAGS)
+LDFLAGS := -L$(lib-prefix)/lib -lblimps $(LDFLAGS)
+
+all: bin
+
+bin: $(BIN)
+
+clean:
+	rm -f *.o $(BIN)
+
+distclean: clean
+	rm -f deps.mk
+
+install:
+	mkdir -p $(DESTDIR)$(prefix)/lib/sift/bin && \
+	cp -a \
+		choose_seqs_via_psiblastseedmedian \
+		clump_output_alignedseq \
+		consensus_to_seq \
+		info_on_seqs \
+		psiblast_res_to_fasta_dbpairwise \
+		seqs_from_psiblast_res \
+		$(DESTDIR)$(prefix)/lib/sift/bin/
+
+help:
+
+.PHONY: all bin clean distclean help install
+
+deps.mk:
+	$(CC) -MM $(CFLAGS) $(addsuffix .c,$(BIN)) > $@
+
+include deps.mk
+
+# vim:ai:
Index: sift/src/Matrix_Info.c
===================================================================
--- sift.orig/src/Matrix_Info.c
+++ sift/src/Matrix_Info.c
@@ -384,7 +384,7 @@ void output_Matrix_Info (Matrix_Info* in
   int pos_observed_tolerant, pos_error_tn, pos_notobserved_intolerant,
 	pos_error_fp, pos_total;
   char c; int l;
-  int print_correct, print_incorrect, print_X, print_type;
+  int print_correct, print_incorrect, print_X, print_type = 0;
 
 	print_correct = 1; print_incorrect = 2; print_X =3; 
   error_fp=0; error_tn = 0, total=0;
@@ -520,7 +520,7 @@ void output_Matrix_Info (Matrix_Info* in
                                                 pos_error_tn++;
 						error_tn++;
 						 print_type = print_incorrect;
-                                        } else {printf ("something is wrong\n");}
+                                        } else {fprintf( stderr, "%s:%d something is wrong\n", __FILE__, __LINE__ );}
                                 	break;
 				case 2: /* +- phenotype is - */
                                         if ( matrix->weights[aa_atob[(int)c]][l] > 0.0) {
@@ -531,7 +531,7 @@ void output_Matrix_Info (Matrix_Info* in
                                                 pos_notobserved_intolerant++;
 						notobserved_intolerant++; 
 						print_type = print_correct;
-                                        }else {printf ("what?\n"); }
+                                        }else {fprintf( stderr, "%s:%d what?\n", __FILE__, __LINE__ ); }
                                         break;
                         	default:
 					fprintf(errorfp, "severity option error %d\n", severity_option);
@@ -612,7 +612,7 @@ void output_Matrix_Info_float (Matrix_In
   int pos_observed_tolerant, pos_error_tn, pos_notobserved_intolerant,
         pos_error_fp, pos_total;
   char c; int l;
-  int print_correct, print_incorrect, print_X, print_type;
+  int print_correct, print_incorrect, print_X, print_type = 0;
 
   print_correct = 1; print_incorrect = 2; print_X =3;
   error_fp=0; error_tn = 0, total=0;
@@ -735,7 +735,7 @@ void output_Matrix_Info_float (Matrix_In
                                                 pos_error_tn++;
                                                 error_tn++;
                                                  print_type = print_incorrect;
-                                        } else {printf ("something is wrong\n");}
+                                        } else {fprintf( stderr, "%s:%d something is wrong\n", __FILE__, __LINE__ );}
                                         break;
                                 case 2: /* +- phenotype is - */
                                         if ( matrix->weights[aa_atob[(int) c]][l] > 0.01) {
@@ -847,10 +847,10 @@ void read_matrix_20_aa_body(mfp, matrix)
    */
 
   len = 0;
-printf ("matrix width %d\n", matrix->width);
+fprintf( stderr, "matrix width %d\n", matrix->width);
   /* skip the letters */
   fgets(Buffer, LARGE_BUFF_LENGTH, mfp);
-printf ("skiipping letters %s\n", Buffer);
+fprintf( stderr, "skipping letters %s\n", Buffer);
   /* read in the numbers */
   while (fgets(Buffer, LARGE_BUFF_LENGTH, mfp) &&
          !blank_line(Buffer) &&
@@ -923,7 +923,7 @@ int
 min_aa_in_column (Matrix * matrix, int pos)
 {
 	double min;
-	int min_aa, aa;
+	int min_aa = 0, aa;
 
 	min = 1000;
 	for (aa = 0; aa < AAS; aa++) {
@@ -1038,7 +1038,7 @@ void output_Matrix_Info_float_uncertaint
   int pos_observed_tolerant, pos_error_tn, pos_notobserved_intolerant,
         pos_error_fp, pos_total;
   char c; int l;
-  int print_correct, print_incorrect, print_X, print_type;
+  int print_correct, print_incorrect, print_X, print_type = 0;
   int uncertain_aa, uncertain_tolerant, uncertain_intolerant;
   int pos_uncertain_aa, pos_uncertain_tolerant, pos_uncertain_intolerant;
   int total_correct, total_counted;
Index: sift/src/PN_blocks.c
===================================================================
--- sift.orig/src/PN_blocks.c
+++ sift/src/PN_blocks.c
@@ -668,7 +668,7 @@ remove_seq0_Xes_from_block (Block* block
         Block* newblock;
         Residue* residue_pointer;
         int i, pos;
-	int non_X_length, newblockpos;
+	int non_X_length, newblockpos = 0;
 
         non_X_length = 0;
 	for (i=0; i < block->width; i++) {
Index: sift/src/PN_convert.c
===================================================================
--- sift.orig/src/PN_convert.c
+++ sift/src/PN_convert.c
@@ -1,6 +1,6 @@
 /* (C) Copyright 2000, Fred Hutchinson Cancer Research Center */
 /* Use, modification or distribution of these programs is subject to */
-/* the terms of the non-commercial licensing agreement in license.h. */   
+/* the terms of the non-commercial licensing agreement in license.h. */
 
 /* convert.c: functions for different methods of converting a block into a */
 /*            matrix */
@@ -8,8 +8,8 @@
 /* Change log information is at the end of the file. */
 /* PN_convert.c copied from convert.c to allow for large block widths */
 
-/*07-13-00 changed in counts() so that if an X is not observed, counts is NOT 
-evenly distributed throughout.  For gaps, counts is evenly distributed 
+/*07-13-00 changed in counts() so that if an X is not observed, counts is NOT
+evenly distributed throughout.  For gaps, counts is evenly distributed
 throughout.
 */
 /*	system headers not in global.h */
@@ -18,6 +18,9 @@ throughout.
 
 #include <assert.h>
 #include <math.h>
+#include <errno.h>
+#include <string.h>
+#include "sift_blimps.h"
 #include "PN_convert.h"
 #include "queue.c"
 #include "SortList.c"
@@ -45,7 +48,7 @@ struct float_qij *Qij;
 void pb_weights(/*block*/);
 
 /*
- * Matrix construction methods.  Build matricies from blocks and 
+ * Matrix construction methods.  Build matricies from blocks and
  * sequence weights.
  */
 
@@ -131,7 +134,7 @@ Matrix *PN_block_to_matrix(block, conver
     break;
   case 22:  /* instead of adding RTot * diffaas*/
 	    /* as in option 20, add RTot * (diffass - 1), thus not adding*/
-	    /* any pseudocounts to completely conserved columns */ 
+	    /* any pseudocounts to completely conserved columns */
 	    /* qij's as pseudocounts */
     PN_altschul_data_dependent_conversion_method (block, matrix, 22);
     break;
@@ -139,13 +142,13 @@ Matrix *PN_block_to_matrix(block, conver
 	   /* exp. pseudocounts*/
     PN_altschul_data_dependent_conversion_method (block, matrix, 23);
     break;
-  case 24: 
+  case 24:
     PN_altschul_data_dependent_conversion_method (block, matrix, 24);
     break;
   case 25:
     PN_altschul_data_dependent_conversion_method (block, matrix, 25);
     break;
-  case 26: 
+  case 26:
     PN_altschul_data_dependent_conversion_method (block, matrix, 26);
     break;
   case 27:
@@ -170,7 +173,7 @@ Matrix *PN_block_to_matrix(block, conver
 	qij_matrix_profile (block, matrix, Qij);
 	break;
    default: /* the default case */
-    sprintf(ErrorBuffer, 
+    sprintf(ErrorBuffer,
 	    "Invalid block to matrix conversion method specified, %d.",
 	    conversion_method);
     ErrorReport(WARNING_ERR_LVL);
@@ -211,10 +214,10 @@ Matrix *PN_block_to_matrix_RTot_par(bloc
 	    ErrorReport(WARNING_ERR_LVL);
 
 	 sprintf(ErrorBuffer,     /* ^^^^----------------vvvvvvvvvvvvvvvvvvv */                    "Using the default conversion method of Altschul's data-dependent method.\n");
- 
+
 	 ErrorReport(WARNING_ERR_LVL);
     	altschul_data_dependent_conversion_method(block, matrix, 0);
-  } 
+  }
 
   /* return the matrix */
   return matrix;
@@ -239,11 +242,11 @@ static struct working *make_col()
 {
   struct working *col;
   int aa;
-  
+
   CheckMem(
 	   col = (struct working *) malloc(sizeof(struct working))
 	   );
-  
+
   col->totcnt = col->totreg = 0.0;
   for (aa=0; aa < AASALL; aa++) {
     col->cnt[aa] = col->reg[aa] = 0.0;
@@ -260,7 +263,7 @@ struct work_pssm *PN_make_pssm(int lengt
   struct work_pssm *pssm;
   int pos, aa;
   double* double_pointer;
- 
+
   CheckMem(
 	   pssm = (struct work_pssm *) malloc(sizeof(struct work_pssm))
 	   );
@@ -268,7 +271,7 @@ struct work_pssm *PN_make_pssm(int lengt
 	double_pointer = (double*) calloc (length * AASALL, sizeof (double));
 	pssm->value = (double **) calloc (length, sizeof (double*));
 	pssm->sum = (double *) calloc (length, sizeof (double));
- 
+
   for (pos = 0; pos < length; pos++) {
 	pssm->value[pos] = &(double_pointer[pos * AASALL]);
   }
@@ -304,7 +307,7 @@ void counts_no_gaps (Block* block, struc
     }
   } /* end of for */
 } /* end of subroutine */
-	
+
 /*======================================================================*/
 static void counts(block, col, pos)
      Block *block;
@@ -314,7 +317,7 @@ static void counts(block, col, pos)
   int seq, aa, aa1;
 
   col->totcnt = col->totraw = col->totreg = 0.0;
-  for (aa = 0; aa < AASALL; aa++) { 
+  for (aa = 0; aa < AASALL; aa++) {
     col->cnt[aa] = col->raw[aa] = col->reg[aa] = 0.0;
   }
 
@@ -356,7 +359,7 @@ static int count_residues(col)
      struct working *col;
 {
   int aa, nr;
-  
+
   nr = 0;
   for (aa = 1; aa < AAS; aa++) {
     if (col->cnt[aa] > 0.0) nr++;
@@ -373,9 +376,9 @@ static void pseudo_alts(col, qij, beta)
      double beta;
 {
   int aa, row;
- 
+
    col->totreg = 0.0;
- 
+
   /*---------- get the pseudo counts -------------------------------*/
   for (aa=1; aa < AAS; aa++) {
     col->reg[aa] = 0.0;
@@ -393,7 +396,7 @@ static void pseudo_alts(col, qij, beta)
 int
 find_min_aa_in_pssm (struct working* col, struct work_pssm* pssm, const int pos)
 {
-	int aa, min_aa;
+	int aa, min_aa = 0;
 	double min_value;
 
 	min_value = 1000;
@@ -420,7 +423,7 @@ void range_on_pssm (struct working* col,
 	min_aa = find_min_aa_in_pssm(col, pssm, pos);
         min = pssm->value[pos][min_aa] * threshold;
         for (aa = 0; aa < AAS; aa++) {
-                pssm->value[pos][aa] -= min; 
+                pssm->value[pos][aa] -= min;
 	}
 
 
@@ -436,14 +439,14 @@ scoring_matrix_profile (Block* block, Ma
 	FILE* matrixfp;
 	char matrix_file[LARGE_BUFF_LENGTH];
 	int pos, original_aa, aa;
-	char* blimps_dir;
+	const char* blimps_dir;
 
 	/* read matrix */
-        blimps_dir = getenv ("BLIMPS_DIR");
+        blimps_dir = getblimpsdir();
 	sprintf (matrix_file, "%s/docs/blosum62.bla.new", blimps_dir);
 
-/*	strcpy (matrix_file, "/tuna/blocks_5.0/blosum/blosum62.bla.new");*/ 
- 
+/*	strcpy (matrix_file, "/tuna/blocks_5.0/blosum/blosum62.bla.new");*/
+
    if ( (matrixfp = fopen (matrix_file, "r")) == NULL) {
         fprintf (errorfp, "can't read the matrix %s\n", matrix_file);
         exit (-1);
@@ -502,17 +505,17 @@ static void compute_BZX(frequency, matri
   /*
    * find the partitions of D, N, E, and Q for B and Z
    */
-  part_D = frequency[aa_atob['D']] / 
+  part_D = frequency[aa_atob['D']] /
     ( frequency[aa_atob['D']] + frequency[aa_atob['N']] );
-  part_N = frequency[aa_atob['N']] / 
+  part_N = frequency[aa_atob['N']] /
     ( frequency[aa_atob['D']] + frequency[aa_atob['N']] );
-  part_E = frequency[aa_atob['E']] / 
+  part_E = frequency[aa_atob['E']] /
     ( frequency[aa_atob['E']] + frequency[aa_atob['Q']] );
-  part_Q = frequency[aa_atob['Q']] / 
+  part_Q = frequency[aa_atob['Q']] /
     ( frequency[aa_atob['E']] + frequency[aa_atob['Q']] );
 
     /* fill in the matrix for B, Z, X, gap, stop and non */
-    matrix->weights[aa_atob['B']][col] = 
+    matrix->weights[aa_atob['B']][col] =
       (part_D * matrix->weights[aa_atob['D']][col] +
        part_N * matrix->weights[aa_atob['N']][col]);
     matrix->weights[aa_atob['Z']][col] =
@@ -538,7 +541,7 @@ static void compute_BZX(frequency, matri
 
 /*=========================================================================
       Adds negative minval to give all positive matrix,
-      then multiplies by 99/maxval to give scores ranging from 0 to 99 
+      then multiplies by 99/maxval to give scores ranging from 0 to 99
       NOTE: Not 0 to 100 because "output_matrix" routine might not leave
 	    enough space.
 ===========================================================================*/
@@ -558,7 +561,7 @@ static void positive_matrix(freqs, pssm,
       if (pssm->value[pos][aa] > maxval) maxval = pssm->value[pos][aa];
     }
   }
-   
+
   if (minval < 0.0) {
     factor = 99.0 / (maxval - minval);
   }
@@ -587,7 +590,7 @@ initialize_cutoff_freqs (double* newfreq
 {
 	int aa;
 
-/* old version. changed July 2, 2001 to make it linux compatible 
+/* old version. changed July 2, 2001 to make it linux compatible
 	assert (oldfreqs[MATRIX_AA_WIDTH] != NULL);
 	assert (newfreqs[MATRIX_AA_WIDTH] != NULL);
 */
@@ -602,7 +605,7 @@ initialize_cutoff_freqs (double* newfreq
 
 
 /* pseudocounts from Altschul's 1997 NAR gapped blast and PSI-BLAST paper */
-void 
+void
 psiblast_alts (Block* block, Matrix* matrix, double* freqs, struct float_qij* qij)
 {
 	double beta;
@@ -652,14 +655,14 @@ psiblast_alts (Block* block, Matrix* mat
   free (pssm->value); free(pssm->sum);
   free(pssm);
 
-} 
+}
 
-int 
+int
 find_max_aa_in_col (struct working* col)
 {
-	int aa, max_aa;
+	int aa, max_aa = 0;
 	double max = 0.0;
-	
+
 	for (aa = 1; aa < AAS; aa++) {
 		if (col->cnt[aa] > max) {
 			max = col->cnt[aa];
@@ -673,7 +676,7 @@ find_max_aa_in_col (struct working* col)
 int
 find_max_aa_in_pssm (struct work_pssm* pssm , int pos)
 {
-	int aa, max_aa;
+	int aa, max_aa = 0;
 	double max;
 
 	max = 0.0;
@@ -685,7 +688,7 @@ find_max_aa_in_pssm (struct work_pssm* p
 	}
 	return max_aa;
 }
-		
+
 /*==========================================================================
      Uses Altschul's method of getting pseudo-counts with a qij matrix,
 ===========================================================================*/
@@ -705,12 +708,12 @@ static void PN_make_alts(block, matrix,
   int original_aa;
   struct dirichlet* diric;
   int min_aa, max_aa;
-  double min_freq, no_of_gaps;
+  double min_freq = 0, no_of_gaps;
   double cutoff_freq[MATRIX_AA_WIDTH];
   int rank_matrix[AAS][AAS];
   double number_of_observed_aas;
   int diri_pseudo;
-  char* blimps_dir;
+  const char* blimps_dir;
   char diri_file[LARGE_BUFF_LENGTH];
 
 	if (scale == 30) {
@@ -723,15 +726,15 @@ static void PN_make_alts(block, matrix,
   col = make_col();
   pssm = PN_make_pssm(block->width);
 
-  blimps_dir = getenv ("BLIMPS_DIR");
-  printf ("blimps dir is %s\n", blimps_dir);
-  sprintf (diri_file , "%s/docs/default.diri", blimps_dir);  
-  printf  ("diri file is %s\n", diri_file);
+  blimps_dir = getblimpsdir();
+  fprintf( stderr, "blimps dir is %s\n", blimps_dir);
+  sprintf (diri_file , "%s/docs/default.diri", blimps_dir);
+  fprintf( stderr, "diri file is %s\n", diri_file);
 
 /*  diric = load_diri_file ("/howard2/pauline/src/merge-opt.13compnew"); */
    diric = load_diri_file (diri_file);
   construct_rank_matrix ( rank_matrix);
-printf ("diri file loaded\n");
+fprintf( stderr, "diri file loaded\n");
   /*--------------  Do one position at a time -------------------*/
   for (pos = 0; pos < block->width; pos++)
   {
@@ -745,8 +748,8 @@ printf ("diri file loaded\n");
 	if (itemp ==1 ) {epsilon = 0; }
 	else {
 		original_aa = block->residues[0][pos];
-                dtemp = similarity_dependent_scale_0 ( col, rank_matrix, original_aa);  
-		epsilon = exp (dtemp); 
+                dtemp = similarity_dependent_scale_0 ( col, rank_matrix, original_aa);
+		epsilon = exp (dtemp);
        }
 
     } else if (scale == 20) {
@@ -762,7 +765,7 @@ printf ("diri file loaded\n");
     /*---------- get the pseudo counts -------------------------------*/
     if (scale == 27 || scale == 29 || scale == 30 || scale == 26) {
                 diri_pseudo = TRUE;
-		pseudo_diric (col, diric, epsilon ); 
+		pseudo_diric (col, diric, epsilon );
     } else {
         pseudo_alts(col, qij, epsilon);
     }
@@ -781,8 +784,8 @@ printf ("diri file loaded\n");
 
        /*    Count+pseudo proportions; returned if scale == 20   */
        if (scale <= 20 || scale == 22 || scale == 25 || scale == 23
-                || scale == 26) 
-                               
+                || scale == 26)
+
        {
           pssm->value[pos][aa] = col->cnt[aa] + col->reg[aa];
         if ( (col->totcnt + col->totreg) > 0.0)
@@ -839,7 +842,7 @@ printf ("diri file loaded\n");
   /*       by subtracting the min. value for all positions ----- */
   if (!scale) positive_matrix(freqs, pssm, matrix);
 
-printf ("left make _alts\n");
+fprintf( stderr, "left make _alts\n");
   free (diric);
   free(col);
   free_struct_work_pssm (pssm );
@@ -852,7 +855,7 @@ printf ("left make _alts\n");
 option:  diri_pseudocounts
 	  TRUE: use 13-component Dirichlet
 	  FALSE use BLOSUM62 qij's
-option:  gap: TRUE -allow everything 
+option:  gap: TRUE -allow everything
 	      FALSE - just look at 20 amino acids
 option:  exp (m)
 	 m=0 # of amino acids at pos (default)
@@ -868,9 +871,9 @@ SIFT_prediction (Block* block, int diri_
 	Matrix* pssm;
 
 	pssm =  copy_block_info_to_matrix (block);
-	
+
         if (Qij == NULL) {
-                printf ("qij matrix missng\n");
+                fprintf( stderr, "qij matrix missing\n");
                 exit (-1);
         }
         normalize  (block);
@@ -912,13 +915,13 @@ void calculate_counts (Block* block)
 		for (aa = 1; aa < AAS; aa++) {
 			printf ("pos %d aa %c counts %.2f weight %.2f\n",
 			pos, aa_btoa[aa], col->cnt[aa], col->totcnt);
-		} 
+		}
 	}
 }
 
-static void SIFT_alts(Block* block, Matrix* matrix, double* freqs, 
-		struct float_qij* qij, 
-		int diri_pseudocounts, int gap_option, 
+static void SIFT_alts(Block* block, Matrix* matrix, double* freqs,
+		struct float_qij* qij,
+		int diri_pseudocounts, int gap_option,
 		int exp_option, int subtract_threshold)
 {
   double dtemp, dtemp2, epsilon;
@@ -933,7 +936,7 @@ static void SIFT_alts(Block* block, Matr
   int div_by_max;
   int rank_matrix[AAS][AAS];
   FILE* rfp;
-  char* blimps_dir;
+  const char* blimps_dir;
   char diri_file[LARGE_BUFF_LENGTH];
 
   div_by_max = TRUE;
@@ -941,25 +944,25 @@ static void SIFT_alts(Block* block, Matr
   col = make_col();
   pssm = PN_make_pssm(block->width);
 
-   blimps_dir = getenv ("BLIMPS_DIR");
+   blimps_dir = getblimpsdir();
   sprintf (diri_file , "%s/docs/default.diri", blimps_dir);
 
    diric = load_diri_file (diri_file);
 
   if (exp_option == 1) {
-          
-	construct_rank_matrix ( rank_matrix); 
+
+	construct_rank_matrix ( rank_matrix);
    }
 
   /*--------------  Do one position at a time -------------------*/
-  for (pos = 0; pos < block->width; pos++) 
+  for (pos = 0; pos < block->width; pos++)
   {
     /*-------- count the number of each aa in this position ------------*/
     if (gap_option) {
 	counts (block,col, pos);
     } else {
-	counts_no_gaps(block, col, pos); 
-   } 
+	counts_no_gaps(block, col, pos);
+   }
 
     /*-------- determine total number of pseudo-counts in column ------*/
     epsilon = 0;
@@ -969,24 +972,24 @@ static void SIFT_alts(Block* block, Matr
     if (exp_option == 1) {
 	/*printf ("pos %d :  ", pos); */
 	dtemp = similarity_dependent_scale_0(col, rank_matrix, original_aa);
-	if (itemp == 1) { 
-		epsilon = 0; 
-	} else { 
+	if (itemp == 1) {
+		epsilon = 0;
+	} else {
 		epsilon = exp (dtemp);
 	}
-    } else { 
+    } else {
 	if (itemp == 1) { epsilon = 0; }
         else if (itemp > 7) { epsilon = 1000; }
         else { epsilon = exp( (double) itemp) ; }
- 
-   } 
-	
+
+   }
+
     /*---------- get the pseudo counts -------------------------------*/
     if (diri_pseudocounts) {
-	pseudo_diric (col, diric, epsilon ); 
+	pseudo_diric (col, diric, epsilon );
     } else {
 	pseudo_alts(col, qij, epsilon);
-    } 
+    }
 
     /*---------   Fill in the matrix entries --------------------*/
     pssm->sum[pos] = 0.0;
@@ -1007,11 +1010,11 @@ static void SIFT_alts(Block* block, Matr
        for (aa = 1; aa < AAS; aa++) {
 		pssm->value[pos][aa] /= min_freq;
 	}
-   } 
+   }
 
    if (subtract_threshold) {
 	   for (aa = 1; aa < AAS; aa++) {
-	        pssm->value[pos][aa] -= TOLERANCE_PROB_THRESHOLD; 
+	        pssm->value[pos][aa] -= TOLERANCE_PROB_THRESHOLD;
    	}
    }
 
@@ -1027,7 +1030,7 @@ static void SIFT_alts(Block* block, Matr
 	if (  (matrix->weights[original_aa][pos] < TOLERANCE_PROB_THRESHOLD
 			 && (!subtract_threshold))
 		||
-	      (matrix->weights[original_aa][pos] < 0.0 && subtract_threshold)		
+	      (matrix->weights[original_aa][pos] < 0.0 && subtract_threshold)
 	) {
 		printf ("WARNING!!! Amino acid %c at pos %d in your original sequence had score %.3f and was not allowed by the prediction.<BR>\n", aa_btoa[original_aa], pos + 1, matrix->weights[original_aa][pos]);
 	}
@@ -1060,7 +1063,7 @@ void SIFT_alts_test(Block* block, Matrix
   int div_by_max;
   int rank_matrix[AAS][AAS];
   FILE* rfp;
-  char* blimps_dir;
+  const char* blimps_dir;
   char diri_file[LARGE_BUFF_LENGTH];
   char matrix_file[LARGE_BUFF_LENGTH];
 
@@ -1070,7 +1073,7 @@ void SIFT_alts_test(Block* block, Matrix
   col = make_col();
   pssm = PN_make_pssm(block->width);
   diri_pseudocounts = TRUE;
-   blimps_dir = getenv ("BLIMPS_DIR");
+   blimps_dir = getblimpsdir();
 
    sprintf (diri_file , "%s/docs/default.diri", blimps_dir);
 
@@ -1080,8 +1083,8 @@ void SIFT_alts_test(Block* block, Matrix
 
         construct_rank_matrix ( rank_matrix);
    }
-  sprintf (matrix_file, "%s/docs/default.qij");
-  init_frq_qij_for_matrix (matrix_file); 
+  sprintf (matrix_file, "%s/docs/default.qij", blimps_dir );
+  init_frq_qij_for_matrix (matrix_file);
 
   for (pos = 0; pos < block->width; pos++)
   {
@@ -1200,7 +1203,7 @@ gap_pseudocounts ( struct working* col,
 
 }
 
-void 
+void
 ratio_mutated_to_normal (struct working* col, int standard_aa)
 {
 	int aa;
@@ -1218,7 +1221,7 @@ assign_gap_positive_scores (Matrix* matr
 {
 	int gaps;
 	int seq;
-	
+
 	assert (matrix->block != NULL);
 /* 	printf ("entering gap_positive scores\n"); */
 	gaps = 0;
@@ -1230,7 +1233,7 @@ assign_gap_positive_scores (Matrix* matr
 	if (gaps == 0) { /* printf ("no gaps\n"); */return; }
 	else { /* printf ("gaps %d\n", gaps) ; */
 			assign_positive_scores_to_column (matrix, pos, gaps); }
-	
+
 } /* end of assign_gap_positive_scores */
 
 void
@@ -1242,7 +1245,7 @@ assign_positive_scores_to_column (Matrix
 	for (aa = 1; aa < AAS; aa++) {
 		matrix->weights[aa][pos] = 0.05;
 	}
-} /* end of positive_column */  		 
+} /* end of positive_column */
 
 /*=========================================================================
     Normalize sequence weights to add up to the number of sequences
@@ -1274,11 +1277,11 @@ PN_altschul_data_dependent_conversion_me
     }
 
   /* check to see if the block has sequence weights */
-printf ("enter normalize block\n");
+fprintf( stderr, "enter normalize block\n");
   normalize(block);   /*  Make weights sum to number of sequences */
-printf ("enter PN_make_alts\n");
-  PN_make_alts(block, matrix, frequency, Qij, RTot, scale); 
-printf ("exit PN_make_alts\n");
+fprintf( stderr, "enter PN_make_alts\n");
+  PN_make_alts(block, matrix, frequency, Qij, RTot, scale);
+fprintf( stderr, "exit PN_make_alts\n");
 }
 
 /*
@@ -1304,7 +1307,7 @@ printf ("exit PN_make_alts\n");
  */
 
 
-/* 
+/*
  * basic_matrix_construction
  *   This is the most general matrix construction method.  Each
  *   occurence of a residue in a sequence contributes the sequence
@@ -1316,7 +1319,7 @@ printf ("exit PN_make_alts\n");
  *     X, '-', '*', & non-code scores are read straight from the frequencies
  *     the frequencies for B and Z are ignored, when a B or Z is encountered
  *       it is partitioned between D & N or E & Q.
- *     the matrix scores for B and Z are computed from the matrix scores of 
+ *     the matrix scores for B and Z are computed from the matrix scores of
  *       D & N and E & Q.
  *   Parameters:
  *     Block *block:       the block to be converted
@@ -1360,13 +1363,13 @@ static void basic_matrix_construction(bl
   /*
    * find the partitions of D, N, E, and Q for B and Z
    */
-  part_D = frequency[aa_atob['D']] / 
+  part_D = frequency[aa_atob['D']] /
     ( frequency[aa_atob['D']] + frequency[aa_atob['N']] );
-  part_N = frequency[aa_atob['N']] / 
+  part_N = frequency[aa_atob['N']] /
     ( frequency[aa_atob['D']] + frequency[aa_atob['N']] );
-  part_E = frequency[aa_atob['E']] / 
+  part_E = frequency[aa_atob['E']] /
     ( frequency[aa_atob['E']] + frequency[aa_atob['Q']] );
-  part_Q = frequency[aa_atob['Q']] / 
+  part_Q = frequency[aa_atob['Q']] /
     ( frequency[aa_atob['E']] + frequency[aa_atob['Q']] );
 
 
@@ -1401,13 +1404,13 @@ static void basic_matrix_construction(bl
       /* otherwise, if it is B, partition B between D and N */
       else if (res == aa_atob['B']) {
 	if (frequency[aa_atob['D']] != 0.0) { /* try to protect div by zero */
-	  count[aa_atob['D']] += 
+	  count[aa_atob['D']] +=
 	    (part_D * seq_weight[seq]) / frequency[aa_atob['D']];
 	  total +=
 	    (part_D * seq_weight[seq]) / frequency[aa_atob['D']];
 	}
 	if (frequency[aa_atob['N']] != 0.0) { /* try to protect div by zero */
-	  count[aa_atob['N']] += 
+	  count[aa_atob['N']] +=
 	    (part_N * seq_weight[seq]) / frequency[aa_atob['N']];
 	  total +=
 	    (part_N * seq_weight[seq]) / frequency[aa_atob['N']];
@@ -1416,13 +1419,13 @@ static void basic_matrix_construction(bl
       /* otherwise, if it is Z, partition Z between E and Q */
       else if (res == aa_atob['Z']) {
 	if (frequency[aa_atob['E']] != 0.0) { /* try to protect div by zero */
-	  count[aa_atob['E']] += 
+	  count[aa_atob['E']] +=
 	    (part_E * seq_weight[seq]) / frequency[aa_atob['E']];
 	  total +=
 	    (part_E * seq_weight[seq]) / frequency[aa_atob['E']];
 	}
 	if (frequency[aa_atob['Q']] != 0.0) { /* try to protect div by zero */
-	  count[aa_atob['Q']] += 
+	  count[aa_atob['Q']] +=
 	    (part_Q * seq_weight[seq]) / frequency[aa_atob['Q']];
 	  total +=
 	    (part_Q * seq_weight[seq]) / frequency[aa_atob['Q']];
@@ -1450,7 +1453,7 @@ static void basic_matrix_construction(bl
     }
 
     /* fill in the matrix for B, Z, X, gap, stop and non */
-    matrix->weights[aa_atob['B']][col] = 
+    matrix->weights[aa_atob['B']][col] =
       (part_D * matrix->weights[aa_atob['D']][col] +
        part_N * matrix->weights[aa_atob['N']][col]);
     matrix->weights[aa_atob['Z']][col] =
@@ -1485,7 +1488,7 @@ get_consensus (Matrix* matrix)
 	int pos;
 
 	seq = (Sequence *) malloc (sizeof (Sequence));
-	seq->sequence = (Residue *) calloc (matrix->width, sizeof (Residue) );	
+	seq->sequence = (Residue *) calloc (matrix->width, sizeof (Residue) );
 	strcpy (seq->name, "CONSENSUS");
 	seq->info[0] = '\0';
 	seq->position = 0;
@@ -1495,7 +1498,7 @@ get_consensus (Matrix* matrix)
 	for (pos = 0; pos < matrix->width; pos++) {
 		res = find_max_aa_in_pos (matrix, pos);
 		seq->sequence[pos] = res;
-		printf ("res %c pos %d\n", aa_btoa[seq->sequence[pos]], pos);  
+		//fprintf( stderr, "res %c pos %d\n", aa_btoa[seq->sequence[pos]], pos);
 	}
 	return seq;
 
@@ -1557,7 +1560,7 @@ cumulative_prob (double col[MATRIX_AA_WI
 {
 	int i;
 	double index;
-	
+
 
 	index = 0;
 
@@ -1569,7 +1572,7 @@ cumulative_prob (double col[MATRIX_AA_WI
 }
 
 
-/* returns a state , given a cumulative probability vector P with 
+/* returns a state , given a cumulative probability vector P with
 MATRIX_AA_WIDTH states ranging from 0 to 1 */
 int
 random_pick (double* P)
@@ -1580,13 +1583,13 @@ random_pick (double* P)
 /*	printf ("random pick: "); */
         random =   (double) rand() /  (double) RAND_MAX * 100 ;
         while (random > P[i] && i < MATRIX_AA_WIDTH) {
-/*                printf ("| %.3f, %c |", P[i], aa_btoa[i]); */ 
+/*                printf ("| %.3f, %c |", P[i], aa_btoa[i]); */
 		i++;
         }
 /*	printf ("picked %c with rand %.3f prob %.3f\n", aa_btoa[i], random, P[i]); */
 	if (i == MATRIX_AA_WIDTH) { /* no aa in this column, return an X */
 		return (aa_atob['X']);
-	} else {	
+	} else {
         	return (i);
 	}
 }
@@ -1772,21 +1775,21 @@ int scores[24][24];
 
 /*=====================================================================*/
 /*   This routine initializes the blimps global variables
-        Qij, RTot and frequency[]      
+        Qij, RTot and frequency[]
 Input Parameters .qij file
 		 .out file (Jorja's format, to get frq)
 */
-void init_frq_qij_for_matrix(char qijname[LARGE_BUFF_LENGTH]) 
+void init_frq_qij_for_matrix(char qijname[LARGE_BUFF_LENGTH])
 {
    FILE *fqij;
- 
+
    if (Qij != NULL) {
 		free (Qij);
 		Qij = NULL;
    }
    if ( (fqij=fopen(qijname, "r")) != NULL)
    { Qij = load_qij(fqij); fclose(fqij);  }
-   else {printf ("could not open %s qij file\n", qijname); }
+   else {fprintf( stderr, "could not open %s qij file\n", qijname); }
    RTot = LOCAL_QIJ_RTOT;
 }  /* end of frq_qij */
 
@@ -1843,7 +1846,7 @@ void pseudo_diric(struct working* col, s
       {
          col->probj[j] = log(diric->q[j]) + col->probn[j] - denom;
 /*       printf("j=%d probn[j]=%f probj[j]=%f\n",
-             j, exp(col->probn[j]), exp(col->probj[j])); */ 
+             j, exp(col->probn[j]), exp(col->probj[j])); */
 
       }
 
@@ -1852,31 +1855,31 @@ void pseudo_diric(struct working* col, s
 	 for (aa = 1; aa < AAS; aa++)
       {
          for (j = 0; j < diric->ncomp; j++) {
-		col->reg[aa] +=   
+		col->reg[aa] +=
 				(exp(col->probj[j]) * diric->alpha[j][aa]);
-			/*  * epsilon); */ 
+			/*  * epsilon); */
             }
 	 col->totreg += col->reg[aa];
 	}
 /*	printf (" total prob. %.2f\n", col->totreg); */
-/* scale dirichlet to probabilities */      
+/* scale dirichlet to probabilities */
 	total = 0.0;
 	for (aa = 1; aa < AAS; aa ++) {
 		col->reg[aa] /= col->totreg;
 		total += col->reg[aa];
 	}
 /*	printf ("%.2f total\n", total); */
- 
+
 /*	printf ("%c %.3f", aa_btoa[aa], col->reg[aa]);
-	if (col->cnt[aa] > 0.0) {printf ("*\n"); } else { printf ("\n");} */ 
+	if (col->cnt[aa] > 0.0) {printf ("*\n"); } else { printf ("\n");} */
 
 }  /* end of pseudo_diric */
 
 /* aa observed in block that has small pseudocount */
-int 
+int
 find_min_aa_in_col_reg (struct working* col)
 {
-        int aa, min_aa;
+        int aa, min_aa = 0;
         double min;
 
         min = 10000;
@@ -1890,10 +1893,10 @@ find_min_aa_in_col_reg (struct working*
         return min_aa;
 }
 
-char 
+char
 find_min_aa_in_col (struct working* col)
 {
-	int aa, min_aa;
+	int aa, min_aa = 0;
 	double min;
 
 	min = 1000;
@@ -1921,7 +1924,7 @@ struct dirichlet *load_diri_file (char f
         struct dirichlet *diric;
 	double denom;
 	double background_frequency[21];
-    printf ("filename is %s\n", filename);
+    fprintf( stderr, "filename is %s\n", filename);
 	if ((fin = fopen (filename, "r")) == NULL) {
 		fprintf (errorfp, "Cannot open dirichlet file %s\n", filename);
 		exit (-1);
@@ -1972,7 +1975,7 @@ struct dirichlet *load_diri_file (char f
 	for (aa = 1; aa < AAS; aa++) {
 		denom = 0.0;
 		for (i = 0; i < numc; i++) {
-			denom += (diric->q[i] * diric->alpha[i][aa] / 
+			denom += (diric->q[i] * diric->alpha[i][aa] /
 							diric->altot[i]);
 		}
 		background_frequency[aa] = denom;
@@ -1980,7 +1983,7 @@ struct dirichlet *load_diri_file (char f
 	for (i =0; i < numc; i++) {
 	/*	printf ("Component %d ratio of aa relative to background \n", i); */
 		for (aa = 1; aa < AAS; aa++) {
-			diric->frequency_to_background_ratio[i][aa] =  
+			diric->frequency_to_background_ratio[i][aa] =
 	diric->alpha[i][aa]/(diric->altot[i] * background_frequency[aa]);
 		}
 	}
@@ -1990,7 +1993,7 @@ struct dirichlet *load_diri_file (char f
 
 
 
-void 
+void
 cutoff_target_freq (double* cutoff_freq, int original_aa, struct float_qij *qij)
 {
 
@@ -2008,8 +2011,9 @@ cutoff_target_freq (double* cutoff_freq,
 void init_frq_qij()
 {
    FILE *fqij;
-   char *blimps_dir, qijname[80], frqname[80];
-   blimps_dir = getenv ("BLIMPS_DIR");
+   const char *blimps_dir;
+   char qijname[80], frqname[80];
+   blimps_dir = getblimpsdir();
 	/*"/howard/servers/blimps/blimps";*/ /* getenv("BLIMPS_DIR"); */
    if (blimps_dir != NULL)
    {
@@ -2025,6 +2029,13 @@ void init_frq_qij()
    Qij = NULL;
    if ( (fqij=fopen(qijname, "r")) != NULL)
    { Qij = load_qij(fqij); fclose(fqij);  }
+   else
+   {
+      fprintf (stderr, "Can't open file '%s': %s.\n", qijname, strerror(errno) );
+      fprintf (stderr, "Current dir:\n", qijname );
+      system( "pwd" );
+      exit (-1);
+   }
    RTot = LOCAL_QIJ_RTOT;
 }  /* end of frq_qij */
 
@@ -2043,7 +2054,7 @@ similarity_dependent_scale_0 (struct wor
   original_aa = max_aa; /* change to max aa 10/24/00 */
 
   n = 0;
-  sum = 0.0; 
+  sum = 0.0;
   for (aa = 1; aa <= 20 ; aa++) {
         if (col->cnt[aa] > 0.0) {
                 n++;
@@ -2065,20 +2076,20 @@ construct_rank_matrix (int rank_matrix[A
 	struct float_qij* rank_Qij;
 	FILE* qijfp;
 	int aa;
-	char* blimps_dir;
+	const char* blimps_dir;
 	char filename[80];
 
-	blimps_dir = getenv ("BLIMPS_DIR");
+	blimps_dir = getblimpsdir();
 	sprintf (filename, "%s/docs/default.sij", blimps_dir);
-	
+
 	if ( (qijfp = fopen (filename, "r"))
 							 == NULL) {
 		fprintf (stderr, "Can't opeen blosum62.sij file in construct_rank_matrix looked in %s\n", blimps_dir);
 		exit (-1);
 	}
- 
+
 	rank_Qij = load_qij (qijfp);
-        fclose (qijfp);	
+        fclose (qijfp);
 
         for (original_aa = 0; original_aa < AAS; original_aa++) {
                 copy_values_to_ranklist (aalist, original_aa, rank_Qij);
@@ -2096,14 +2107,14 @@ assign_ranks (int rank_matrix[AAS][AAS],
 
         for (rank = 0; rank < AAS; rank++) {
                 aa = aalist[rank].aa;
-                rank_matrix[original_aa][aa] = rank + 1; 
+                rank_matrix[original_aa][aa] = rank + 1;
 			/* otherwise rank starts from 0
                              to AAS-1 rather than 1 to AAS */
         }
 }
 
 void
-copy_values_to_ranklist (struct rank_cell aalist[AAS], int original_aa, 
+copy_values_to_ranklist (struct rank_cell aalist[AAS], int original_aa,
 			struct float_qij* rankQij) {
 
         int aa;
@@ -2111,8 +2122,8 @@ copy_values_to_ranklist (struct rank_cel
 	/* 01/03/00 rankQij has -1 for gap, which will mess up
 	ranks.  Assign it a really negative value so it will have
 	the lowest rank */
-	
-	aalist[0].aa = 0; 
+
+	aalist[0].aa = 0;
 	aalist[0].value = -50;
 
         for (aa = 1; aa < AAS; aa++) {
@@ -2180,7 +2191,7 @@ number_of_real_aminoacids (Block* block,
 
 	count = 0;
 	for (seq = 0; seq < block->num_sequences; seq++) {
-		if (block->residues[seq][pos] >= 1 && 
+		if (block->residues[seq][pos] >= 1 &&
 			block->residues[seq][pos] <=20) {
 			count++;
 		}
@@ -2189,7 +2200,7 @@ number_of_real_aminoacids (Block* block,
 } /* end of number_of_real_aminoacids */
 
 double*
-calculate_basic_aa_fraction (Block* block) 
+calculate_basic_aa_fraction (Block* block)
 {
 	double* fraction_stored;
 	double percent_obsaa_div_seq;
@@ -2219,15 +2230,15 @@ calculate_basic_aa_stored (Block* block)
 }
 
 void
-free_struct_work_pssm (struct work_pssm* pssm) 
+free_struct_work_pssm (struct work_pssm* pssm)
 {
 	free (pssm->value[0]); /* holy cow, Ithink this works. This points
                                to space allocated in double pointer*/
 	free (pssm->value); /* this points to an array of block length .
-		  each cell points to values allocated by double pointer*/ 
+		  each cell points to values allocated by double pointer*/
 	free (pssm->sum);
 	free (pssm);
-	pssm = NULL;	
+	pssm = NULL;
 
 }
 
Index: sift/src/Protdist.c
===================================================================
--- sift.orig/src/Protdist.c
+++ sift/src/Protdist.c
@@ -124,7 +124,7 @@ seq_with_max_dist (struct score_pair* pa
 {
 	double max_dist;
 	int i;
-	int index_of_max_seq;
+	int index_of_max_seq = 0;
 
 	max_dist = 0;
 
Index: sift/src/Psiblast.c
===================================================================
--- sift.orig/src/Psiblast.c
+++ sift/src/Psiblast.c
@@ -1,6 +1,6 @@
 /* (C) Copyright 2000, Fred Hutchinson Cancer Research Center */
 /* Use, modification or distribution of these programs is subject to */
-/* the terms of the non-commercial licensing agreement in license.h. */   
+/* the terms of the non-commercial licensing agreement in license.h. */
 
 #ifndef _PSIBLAST_C_
 #define _PSIBLAST_C_
@@ -20,7 +20,7 @@ read_psiblast_entry ( FILE* fp, char Buf
 /* passing first line in entry that begins with >seq_id in Buffer */
         Aligned_Pair* alignment, *new_alignment;
 	Aligned_Pair * cursor;
-	
+
         alignment = initialize_Aligned_Pair (Buffer);
 	cursor = alignment;
 
@@ -35,14 +35,14 @@ read_psiblast_entry ( FILE* fp, char Buf
 		new_alignment =initialize_Aligned_Pair (Buffer);
 		cursor->next = new_alignment;
 		cursor = new_alignment;
-	 	reading_alignment_at_score_line (cursor, Buffer, fp);		
+	 	reading_alignment_at_score_line (cursor, Buffer, fp);
 	}
 	return (alignment);
 
 } /* end of read_psiblast entry*/
 
 Aligned_Pair*
-initialize_Aligned_Pair (char Buffer[LARGE_BUFF_LENGTH]) 
+initialize_Aligned_Pair (char Buffer[LARGE_BUFF_LENGTH])
 {
 	Aligned_Pair* alignment;
 	char* stringptr;
@@ -70,7 +70,7 @@ initialize_Aligned_Pair (char Buffer[LAR
 } /* end of initialize_Aligned_Pair*/
 
 void
-reading_alignment_at_score_line (Aligned_Pair* alignment, 
+reading_alignment_at_score_line (Aligned_Pair* alignment,
 				char Buffer[LARGE_BUFF_LENGTH], FILE* fp)
 {
         char* buff;
@@ -167,7 +167,7 @@ void process_sequence_line(Sequence* seq
   /* Note: The buff may have characters that are not residues, this is */
   /*       ok, there will just be a little extra space */
   while ((saved_length + estimated_length) > seq->max_length) {
-    resize_sequence(seq); 
+    resize_sequence(seq);
   }
 
   /* temporarily increase the room to put the sequence into */
@@ -194,23 +194,23 @@ free_aligned_pair (Aligned_Pair* alignme
 		free (alignment);
 	}
 	alignment = NULL;
-	
+
 }
 
-void block_to_checkpoint_file (Block* block, FILE* fp, FILE* outfp) 
+void block_to_checkpoint_file (Block* block, FILE* fp, FILE* outfp)
 {
 
 	Matrix* pssm;
-	long ltemp;
+	int32_t ltemp;
 	int pos;
 	char ctemp[12];
 	int aa; double dtemp;
-printf ("in block to checkpoint file\n");
+fprintf( stderr, "in block to checkpoint file\n");
 	pssm = PN_block_to_matrix (block, 20);
-/*	print_matrix (pssm);*/ 
-printf ("after conversion\n");
-	ltemp = (long) block->width;
-	fwrite (&ltemp, sizeof (long), 1, fp);
+/*	print_matrix (pssm);*/
+fprintf( stderr, "after conversion\n");
+	ltemp = (int32_t) block->width;
+	fwrite (&ltemp, sizeof(ltemp), 1, fp);
 	for (pos = 0; pos < block->width; pos++)
 	{
 		ctemp[0] = toupper (aa_btoa[block->sequences[0].sequence[pos]]);
@@ -313,7 +313,7 @@ get_top_seq (char next_seq_to_add[], Has
 
 
 	if (read_psiblast_header_until_first_no_error (psiblastfp, FALSE) == -1) {
-		next_seq_to_add[0] = '\0'; 
+		next_seq_to_add[0] = '\0';
 	} else {
 
 	fgets (line, LINE_LEN, psiblastfp);
@@ -322,24 +322,24 @@ get_top_seq (char next_seq_to_add[], Has
 	assert ( strlen(strptr) < KEY_WIDTH);
 	strncpy (name, strptr, KEY_WIDTH);
 	name[KEY_WIDTH-1] = '\0';
-	} 
+	}
 	while (Exists(name, namehash)) {
-		fgets (line, LINE_LEN, psiblastfp); 
+		fgets (line, LINE_LEN, psiblastfp);
 		strptr = strtok(line, " \t\n");
         	if (strptr == NULL) {
-			printf ("entered no line\n");
+			fprintf( stderr, "entered no line\n");
 			name[0] = '\0';
 		} else {
 			strncpy (name, strptr, KEY_WIDTH);
 			/* 	strptr[strlen(strptr)] = '\0' ; */
-			name[KEY_WIDTH-1] = '\0';      
+			name[KEY_WIDTH-1] = '\0';
 		}
 	}
 	strncpy (next_seq_to_add, name, KEY_WIDTH);
 	next_seq_to_add[KEY_WIDTH -1] = '\0';
 
 	fclose (psiblastfp);
- 
+
 } /* end get_top_seq */
 
 char* substring(Sequence* seq, int begin, int len)
Index: sift/src/psiblast_res_to_fasta_dbpairwise.c
===================================================================
--- sift.orig/src/psiblast_res_to_fasta_dbpairwise.c
+++ sift/src/psiblast_res_to_fasta_dbpairwise.c
@@ -32,8 +32,8 @@ Output2: Sequences found with residues f
 #include "Protdist.c"
 #include "Psiblast.c"
 
-#define MAXSEQ 400 /* Maximum number of sequences */
-#define MIN_SEQ 5
+#define MAXSEQ 20000 /* Maximum number of sequences */
+#define MIN_SEQ 2
 #define LINE_LEN 800
 #define FALSE 0
 #define TRUE 1
@@ -113,17 +113,17 @@ int main
 
 	free_seqs (alignedseqs, nseqs);
         if (nseqs < MIN_SEQ) {
-                fprintf (errorfp, "ERROR! Not enough sequences (only %d) found by the PSI-BLAST search! Program terminated.\n", nseqs -1);
-                exit (-1);
+                fprintf (stderr,  "Warning: few sequences (only %d, MIN_SEQ = %d) found by the PSI-BLAST search.\n", nseqs -1, MIN_SEQ);
+                fprintf (errorfp, "Warning: few sequences (only %d, MIN_SEQ = %d) found by the PSI-BLAST search.\n", nseqs -1, MIN_SEQ);
+                exit(0);
         } else if (nseqs > MAXSEQ) {
-                fprintf (errorfp, "WARNING: exceeded maximum # of alignments in processing PSIBLAST result\n");
-                exit (-1);
+                fprintf (stderr,  "ERROR: exceeded maximum number of alignments (%d > %d) in processing PSIBLAST result\n", nseqs, MAXSEQ );
+                fprintf (errorfp, "ERROR: exceeded maximum number of alignments (%d > %d) in processing PSIBLAST result\n", nseqs, MAXSEQ );
+                exit(-1);
         }
 	fclose (errorfp);
 	rm_file (errorfilename);
 	exit (0);
-
-
 } /* end main */
 
 void getargs (int argc, char* argv[], FILE** psiblastfp, 
@@ -201,6 +201,7 @@ read_sequence_from_filename (char filena
 	Sequence* seq;
 
 	if ((fp = fopen (filename, "r")) == NULL) {
+		fprintf (stderr,  "couldn't open %s in processing PSIBLAST results\n", filename);
 		fprintf (errorfp, "couldn't open %s in processing PSIBLAST results\n", filename);
 		exit (-1);
 	}
Index: sift/src/seqs_from_psiblast_res.c
===================================================================
--- sift.orig/src/seqs_from_psiblast_res.c
+++ sift/src/seqs_from_psiblast_res.c
@@ -23,7 +23,7 @@ and prints out the complete sequence fro
 #include "stringhash.c"
 #include "Psiblast.c"
 
-#define MAXSEQ 400 /* Maximum number of sequences */
+#define MAXSEQ 20000 /* Maximum number of sequences */
 #define LINE_LEN 800
 #define FALSE 0
 #define TRUE 1
Index: sift/src/sift_blimps.h
===================================================================
--- /dev/null
+++ sift/src/sift_blimps.h
@@ -0,0 +1,12 @@
+#ifndef SIFT_BLIMPS_H
+#define SIFT_BLIMPS_H
+
+inline const char*		getblimpsdir()
+{
+	const char* blimps_dir = getenv ("BLIMPS_DIR");
+	if( blimps_dir == NULL ) blimps_dir = __MAKE_PREFIX__ "/share/sift/blimps";
+
+	return blimps_dir;
+}
+
+#endif
Index: sift/src/SortList.c
===================================================================
--- sift.orig/src/SortList.c
+++ sift/src/SortList.c
@@ -13,7 +13,7 @@ int
 lowest_scoring_aa (Matrix* matrix, int pos)
 {
 	struct working* col;
-	int aa, min_aa;
+	int aa, min_aa = 0;
 	double min;
 
 	min = 10000;
Index: sift/SIFT_for_submitting_fasta_seq.csh.pod
===================================================================
--- /dev/null
+++ sift/SIFT_for_submitting_fasta_seq.csh.pod
@@ -0,0 +1,55 @@
+=head1 NAME
+
+SIFT_for_submitting_fasta_seq.csh - predict effect of an amino acid substitution on protein function
+
+=head1 SYNOPSIS
+
+SIFT_for_submitting_fasta_seq.csh <FASTA_FILE> <BLAST_DB> <SUBSTITUTIONS_FILE> [BLAST_PROCESSORS]
+
+=head1 DESCRIPTION
+
+SIFT predicts whether an amino acid substitution affects protein function based on sequence homology and the physical properties of amino acids.
+
+Results are stored in F<< ./<seq_file>.SIFTprediction >>.
+
+This program is used for FASTA input and is part of the SIFT suite.
+
+=head1 OPTIONS
+
+=over
+
+=item <FASTA_FILE>
+
+Protein sequence in fasta format.
+
+=item <BLAST_DB>
+
+Protein database to search.  These sequences are assumed to be functional.
+
+=item <SUBSTITUTIONS_FILE>
+
+File of substitutions to be predicted.  See /usr/share/doc/sift/examples/lacI.subst for an example of the format.  If you give '-', scores for all mutations of the entire protein sequence are printed.
+
+=item [BLAST_PROCESSORS]
+
+Number of processors/cores to use when running blast via the B<-a> argument.
+
+=back
+
+=head1 EXAMPLES
+
+ SIFT_for_submitting_fasta_seq.csh /usr/share/doc/sift/examples/lacI.fasta [BLAST_DB] /usr/share/doc/sift/examples/lacI.subst
+
+=head1 AUTHOR
+
+Pauline Ng
+
+=head1 COPYRIGHT AND LICENSE
+
+(C) Copyright 1993-2001, Fred Hutchinson Cancer Research Center
+
+Noncommercial license. See /usr/share/doc/sift/copyright for details.
+
+=head1 SEE ALSO
+
+L<info_on_seqs>
Index: sift/src/Clumping.c
===================================================================
--- sift.orig/src/Clumping.c
+++ sift/src/Clumping.c
@@ -175,7 +175,7 @@ printf ("entered if statment\n");
  printf ("Clumped %d sequences into %d clusters\n", nseq, iclus);
  /* keep this print because get # of clustersin awk statement */
   free(pairs);
-printf ("exiting clump\n");
+fprintf( stderr, "%s:%d exiting clump\n", __FILE__, __LINE__);
  return clusters_of_seqs;
 }  /* end of cluster */
 
@@ -432,7 +432,7 @@ printf ("entered if statment\n");
  printf ("Clumped %d sequences into %d clusters\n", nseq, i2);
  /* keep this print because get # of clustersin awk statement */
   free(pairs);
-printf ("exiting clump\n");
+fprintf( stderr, "%s:%d exiting clump\n", __FILE__, __LINE__);
  return no_of_clumps;
 }  /* end of cluster */
 
Index: sift/info_on_seqs.pod
===================================================================
--- /dev/null
+++ sift/info_on_seqs.pod
@@ -0,0 +1,46 @@
+=head1 NAME
+
+info_on_seqs - perform final step of SIFT prediction
+
+=head1 SYNOPSIS
+
+info_on_seqs <FASTA_FILE> <SUBSTITUTIONS_FILE> <OUTPUT>
+
+=head1 DESCRIPTION
+
+You should use L<SIFT_for_submitting_fasta_seq.csh> instead of this program unless you know what you are doing.
+
+This program prepares the final SIFT prediction.  In addition to the command line arguments, it expects certain files present I<next to> the FASTA_FILE.
+
+=head1 OPTIONS
+
+=over
+
+=item <FASTA_FILE>
+
+Protein sequence in fasta format.
+
+=item <SUBSTITUTIONS_FILE>
+
+File of substitutions to be predicted.  See /usr/share/doc/sift/examples/lacI.subst for an example of the format.  If you give '-', scores for all mutations of the entire protein sequence are printed.
+
+=item <OUTPUT>
+
+Output file.
+
+=back
+
+=head1 AUTHOR
+
+Pauline Ng
+
+=head1 COPYRIGHT AND LICENSE
+
+(C) Copyright 1993-2001, Fred Hutchinson Cancer Research Center
+
+Noncommercial license. See /usr/share/doc/sift/copyright for details.
+
+=head1 SEE ALSO
+
+L<SIFT_for_submitting_fasta_seq.csh>
+