blast_options.c

ncbi源码

/*
 * ===========================================================================
 * PRODUCTION $Log: blast_options.c,v $
 * PRODUCTION Revision 1000.6  2004/06/01 18:07:31  gouriano
 * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.111
 * PRODUCTION
 * ===========================================================================
 */

/* $Id: blast_options.c,v 1000.6 2004/06/01 18:07:31 gouriano Exp $
 * ===========================================================================
 *
 *                            PUBLIC DOMAIN NOTICE
 *               National Center for Biotechnology Information
 *
 *  This software/database is a "United States Government Work" under the
 *  terms of the United States Copyright Act.  It was written as part of
 *  the author's offical duties as a United States Government employee and
 *  thus cannot be copyrighted.  This software/database is freely available
 *  to the public for use. The National Library of Medicine and the U.S.
 *  Government have not placed any restriction on its use or reproduction.
 *
 *  Although all reasonable efforts have been taken to ensure the accuracy
 *  and reliability of the software and data, the NLM and the U.S.
 *  Government do not and cannot warrant the performance or results that
 *  may be obtained by using this software or data. The NLM and the U.S.
 *  Government disclaim all warranties, express or implied, including
 *  warranties of performance, merchantability or fitness for any particular
 *  purpose.
 *
 *  Please cite the author in any work or product based on this material.
 *
 * ===========================================================================
 */

/** @file blast_options.c
 * @todo FIXME needs file description & doxygen comments
 */

static char const rcsid[] = 
    "$Id: blast_options.c,v 1000.6 2004/06/01 18:07:31 gouriano Exp $";

#include <algo/blast/core/blast_options.h>
#include <algo/blast/core/blast_filter.h>
#include <algo/blast/core/blast_encoding.h>

QuerySetUpOptions*
BlastQuerySetUpOptionsFree(QuerySetUpOptions* options)

{
   sfree(options->filter_string);

   sfree(options);
   return NULL;
}

Int2
BlastQuerySetUpOptionsNew(QuerySetUpOptions* *options)
{
   *options = (QuerySetUpOptions*) calloc(1, sizeof(QuerySetUpOptions));
   
   if (*options == NULL)
      return 1;
   
   (*options)->genetic_code = BLAST_GENETIC_CODE;

   return 0;
}

Int2 BLAST_FillQuerySetUpOptions(QuerySetUpOptions* options,
        Uint1 program, const char *filter_string, Uint1 strand_option)
{
   if (options == NULL)
      return 1;
   
   if (strand_option && 
       (program == blast_type_blastn || program == blast_type_blastx ||
       program == blast_type_tblastx)) {
      options->strand_option = strand_option;
   }

   /* "L" indicates low-complexity (seg for proteins, 
      dust for nucleotides). */
   if (!filter_string || !strcasecmp(filter_string, "T")) {
      if (program == blast_type_blastn)
         options->filter_string = strdup("D");
      else
         options->filter_string = strdup("S");
   } else {
      options->filter_string = strdup(filter_string); 
   }

   return 0;
}

BlastInitialWordOptions*
BlastInitialWordOptionsFree(BlastInitialWordOptions* options)

{

	sfree(options);

	return NULL;
}


Int2
BlastInitialWordOptionsNew(Uint1 program, 
   BlastInitialWordOptions* *options)
{
   *options = 
      (BlastInitialWordOptions*) calloc(1, sizeof(BlastInitialWordOptions));
   if (*options == NULL)
      return 1;

   if (program != blast_type_blastn) {	/* protein-protein options. */
      (*options)->window_size = BLAST_WINDOW_SIZE_PROT;
      (*options)->x_dropoff = BLAST_UNGAPPED_X_DROPOFF_PROT;
   } else {
      (*options)->window_size = BLAST_WINDOW_SIZE_NUCL;
   }

   return 0;
}

Int2
BLAST_FillInitialWordOptions(BlastInitialWordOptions* options, 
   Uint1 program, Boolean greedy, Int4 window_size, 
   Boolean variable_wordsize, Boolean ag_blast, Boolean mb_lookup,
   double xdrop_ungapped)
{
   if (!options)
      return 1;

   /* Ungapped extension is performed in all cases except when greedy
      gapped extension is used */
   if (program != blast_type_blastn) {
      options->ungapped_extension = TRUE;
      options->x_dropoff = BLAST_UNGAPPED_X_DROPOFF_PROT;
   } else if (!greedy) {
      options->ungapped_extension = TRUE;
      options->x_dropoff = BLAST_UNGAPPED_X_DROPOFF_NUCL;
   } else {
      options->ungapped_extension = FALSE;
   }

   if (window_size != 0)
      options->window_size = window_size;
   if (xdrop_ungapped != 0)
      options->x_dropoff = xdrop_ungapped;

   options->variable_wordsize = variable_wordsize;

   if (ag_blast) {
      options->extension_method = eRightAndLeft;
   } else {
      options->extension_method = eRight;
      if (mb_lookup)
         options->container_type = eMbStacks;
      else
         options->container_type = eDiagArray;
   }

   return 0;
}



BlastInitialWordParameters*
BlastInitialWordParametersFree(BlastInitialWordParameters* parameters)

{
	sfree(parameters);
	return NULL;

}

/** Compute the default cutoff expect value for ungapped extensions
 * @param program The blast program type
 * @return The default per-program expect value
 */
static double GetUngappedCutoff(Uint1 program)
{
   switch(program) {
   case blast_type_blastn:
      return UNGAPPED_CUTOFF_E_BLASTN;
   case blast_type_blastp: 
   case blast_type_rpsblast: 
      return UNGAPPED_CUTOFF_E_BLASTP;
   case blast_type_blastx: 
      return UNGAPPED_CUTOFF_E_BLASTX;
   case blast_type_tblastn:
   case blast_type_rpstblastn:
      return UNGAPPED_CUTOFF_E_TBLASTN;
   case blast_type_tblastx:
      return UNGAPPED_CUTOFF_E_TBLASTX;
   }
   return 0;
}

Int2
BlastInitialWordParametersNew(Uint1 program_number, 
   const BlastInitialWordOptions* word_options, 
   const BlastHitSavingParameters* hit_params, 
   const BlastExtensionParameters* ext_params, BlastScoreBlk* sbp, 
   BlastQueryInfo* query_info, 
   Uint4 subject_length,
   BlastInitialWordParameters* *parameters)
{
   Int2 status = 0;

   /* If parameters pointer is NULL, there is nothing to fill, 
      so don't do anything */
   if (!parameters)
      return 0;

   ASSERT(sbp);
   ASSERT(word_options);

   *parameters = (BlastInitialWordParameters*) 
      calloc(1, sizeof(BlastInitialWordParameters));

   (*parameters)->options = (BlastInitialWordOptions *) word_options;

   ASSERT(sbp->kbp_std[query_info->first_context]);

   (*parameters)->x_dropoff = (Int4)
      ceil(sbp->scale_factor * word_options->x_dropoff * NCBIMATH_LN2/
           sbp->kbp_std[query_info->first_context]->Lambda);

   status = BlastInitialWordParametersUpdate(program_number,
               hit_params, ext_params, sbp, query_info,
               subject_length, *parameters);

   return status;
}

Int2
BlastInitialWordParametersUpdate(Uint1 program_number, 
   const BlastHitSavingParameters* hit_params, 
   const BlastExtensionParameters* ext_params, BlastScoreBlk* sbp, 
   BlastQueryInfo* query_info, Uint4 subj_length,
   BlastInitialWordParameters* parameters)
{
   Blast_KarlinBlk* kbp;
   Boolean gapped_calculation = TRUE;

   ASSERT(sbp);
   ASSERT(hit_params);
   ASSERT(ext_params);
   ASSERT(query_info);

   /* kbp_gap is only non-NULL for gapped searches! */
   if (sbp->kbp_gap) {
      kbp = sbp->kbp_gap[query_info->first_context];
   } else {
      kbp = sbp->kbp_std[query_info->first_context];
      gapped_calculation = FALSE;
   }

   ASSERT(kbp);
   /* For non-blastn programs cutoff score should not be larger than 
      gap trigger. */
   if (gapped_calculation && program_number != blast_type_blastn) {
      parameters->cutoff_score = 
         MIN(ext_params->gap_trigger, hit_params->cutoff_score);
   } else {
      Int4 s2 = 0;
      double e2;
      Int4 qlen;
      /* Calculate score cutoff corresponding to a fixed e-value (0.05);
         If it is smaller, then use this one */
      qlen = query_info->context_offsets[query_info->last_context+1] - 1;
      
      e2 = GetUngappedCutoff(program_number);

      BLAST_Cutoffs(&s2, &e2, kbp, MIN(subj_length, (Uint4) qlen)*subj_length, TRUE, 
                    hit_params->gap_decay_rate);
      s2 *= (Int4)sbp->scale_factor;
      parameters->cutoff_score = MIN(hit_params->cutoff_score, s2);

      if (parameters->x_dropoff != 0 && parameters->cutoff_score != 0) {
         parameters->x_dropoff = 
            MIN(parameters->x_dropoff, parameters->cutoff_score);
      } else if (parameters->cutoff_score != 0) {
         parameters->x_dropoff = parameters->cutoff_score;
      }
   }

   return 0;
}

BlastExtensionOptions*
BlastExtensionOptionsFree(BlastExtensionOptions* options)

{

	sfree(options);

	return NULL;
}



Int2
BlastExtensionOptionsNew(Uint1 program, BlastExtensionOptions* *options)

{
	*options = (BlastExtensionOptions*) 
           calloc(1, sizeof(BlastExtensionOptions));

	if (*options == NULL)
		return 1;

	if (program != blast_type_blastn) /* protein-protein options. */
	{
		(*options)->gap_x_dropoff = BLAST_GAP_X_DROPOFF_PROT;
		(*options)->gap_x_dropoff_final = 
                   BLAST_GAP_X_DROPOFF_FINAL_PROT;
		(*options)->gap_trigger = BLAST_GAP_TRIGGER_PROT;
		(*options)->ePrelimGapExt = eDynProgExt;
	}
	else
	{
		(*options)->gap_trigger = BLAST_GAP_TRIGGER_NUCL;
	}

	return 0;
}

Int2
BLAST_FillExtensionOptions(BlastExtensionOptions* options, 
   Uint1 program, Boolean greedy, double x_dropoff, double x_dropoff_final)
{
   if (!options)
      return 1;

   if (program == blast_type_blastn) {
      if (greedy) {
         options->gap_x_dropoff = BLAST_GAP_X_DROPOFF_GREEDY;
         options->ePrelimGapExt = eGreedyWithTracebackExt;
         options->eTbackExt = eGreedyTbck;
      }	else {
         options->gap_x_dropoff = BLAST_GAP_X_DROPOFF_NUCL;
         options->gap_x_dropoff_final = BLAST_GAP_X_DROPOFF_FINAL_NUCL;
         options->ePrelimGapExt = eDynProgExt;
         options->eTbackExt = eDynProgTbck;
      }
   }

   if (x_dropoff)
      options->gap_x_dropoff = x_dropoff;
   if (x_dropoff_final) {
      options->gap_x_dropoff_final = x_dropoff_final;
   } else {
      /* Final X-dropoff can't be smaller than preliminary X-dropoff */
      options->gap_x_dropoff_final = 
         MAX(options->gap_x_dropoff_final, x_dropoff);
   }

   return 0;

}

Int2 
BlastExtensionOptionsValidate(Uint1 program_number, 
   const BlastExtensionOptions* options, Blast_Message* *blast_msg)

{
	if (options == NULL)
		return 1;

	if (program_number != blast_type_blastn)
	{
		if (options->ePrelimGapExt == eGreedyWithTracebackExt || 
            	options->ePrelimGapExt == eGreedyExt)
		{
			Int4 code=2;
			Int4 subcode=1;
			Blast_MessageWrite(blast_msg, BLAST_SEV_WARNING, code, subcode, 
                            "Greedy extension only supported for BLASTN");
			return (Int2) code;
		}
	}

	return 0;
}

Int2 BlastExtensionParametersNew(Uint1 program_number, 
        const BlastExtensionOptions* options, BlastScoreBlk* sbp,
        BlastQueryInfo* query_info, BlastExtensionParameters* *parameters)
{
   Blast_KarlinBlk* kbp,* kbp_gap;
   BlastExtensionParameters* params;

   /* If parameters pointer is NULL, there is nothing to fill, 
      so don't do anything */
   if (!parameters)
      return 0;

   if (sbp->kbp) {
      kbp = sbp->kbp[query_info->first_context];
   } else {
      /* The Karlin block is not found, can't do any calculations */
      *parameters = NULL;
      return -1;
   }

   if (sbp->kbp_gap) {
      kbp_gap = sbp->kbp_gap[query_info->first_context];
   } else {
      kbp_gap = kbp;
   }
   
   *parameters = params = (BlastExtensionParameters*) 
      malloc(sizeof(BlastExtensionParameters));

   params->options = (BlastExtensionOptions *) options;
   params->gap_x_dropoff = 
      (Int4) (options->gap_x_dropoff*NCBIMATH_LN2 / kbp_gap->Lambda);
   params->gap_x_dropoff_final = (Int4) 
      (options->gap_x_dropoff_final*NCBIMATH_LN2 / kbp_gap->Lambda);

   params->gap_trigger = (Int4)
      ((options->gap_trigger*NCBIMATH_LN2 + kbp->logK) / kbp->Lambda);

   if (sbp->scale_factor > 1.0) {
      params->gap_trigger *= (Int4)sbp->scale_factor;
      params->gap_x_dropoff *= (Int4)sbp->scale_factor;
      params->gap_x_dropoff_final *= (Int4)sbp->scale_factor;
   }
   return 0;
}

BlastExtensionParameters*
BlastExtensionParametersFree(BlastExtensionParameters* parameters)
{
  sfree(parameters);
  return NULL;
}


BlastScoringOptions*
BlastScoringOptionsFree(BlastScoringOptions* options)

{