blast_options.c

ncbi源码

	{
		Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR, code, subcode, 
                         "Non-zero threshold required");
		return (Int2) code;
	}

	if (options->word_size <= 0)
	{
		Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR, code, subcode, 
                         "Word-size must be greater than zero");
		return (Int2) code;
	} else if (program_number == blast_type_blastn && options->word_size < 7)
	{
		Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR, code, subcode, 
                         "Word-size must be 7" 
                         "or greater for nucleotide comparison");
		return (Int2) code;
	} else if (program_number != blast_type_blastn && options->word_size > 5)
	{
		Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR, code, subcode, 
                         "Word-size must be less"
                         "than 6 for protein comparison");
		return (Int2) code;
	}


	if (program_number != blast_type_blastn && 
       options->lut_type == MB_LOOKUP_TABLE)
	{
		Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR, code, subcode, 
                         "Megablast lookup table only supported with blastn");
		return (Int2) code;
	}

   if (options->lut_type == MB_LOOKUP_TABLE && options->word_size < 12 && 
       options->mb_template_length == 0) {
      Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR, code, subcode, 
                         "Word size must be 12 or greater with megablast"
                         " lookup table");
      return (Int2) code;
   }

   if (program_number == blast_type_blastn && 
       options->mb_template_length > 0) {
      if (!DiscWordOptionsValidate(options->word_size,
            options->mb_template_length, options->mb_template_type)) {
         Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR, code, subcode, 
            "Invalid discontiguous template parameters");
         return (Int2) code;
      } else if (options->lut_type != MB_LOOKUP_TABLE) {
         Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR, code, subcode, 
            "Invalid lookup table type for discontiguous Mega BLAST");
         return (Int2) code;
      } else if (options->scan_step != 1 && 
                 options->scan_step != COMPRESSION_RATIO) {
         Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR, code, subcode, 
            "Invalid scanning stride for discontiguous Mega BLAST");
         return (Int2) code;
      }
   }

	return 0;
}

BlastHitSavingOptions*
BlastHitSavingOptionsFree(BlastHitSavingOptions* options)

{
  sfree(options);
  return NULL;
}


Int2 BlastHitSavingOptionsNew(Uint1 program_number, 
        BlastHitSavingOptions* *options)
{
   *options = (BlastHitSavingOptions*) calloc(1, sizeof(BlastHitSavingOptions));
   
   if (*options == NULL)
      return 1;

   (*options)->hitlist_size = 500;
   (*options)->prelim_hitlist_size = 500;
   (*options)->expect_value = BLAST_EXPECT_VALUE;

   /* other stuff?? */
   
   return 0;

}

Int2
BLAST_FillHitSavingOptions(BlastHitSavingOptions* options, 
                           double evalue, Int4 hitlist_size)
{
   if (!options)
      return 1;

   if (hitlist_size)
      options->hitlist_size = options->prelim_hitlist_size = hitlist_size;
   if (evalue)
      options->expect_value = evalue;

   return 0;

}

Int2
BlastHitSavingOptionsValidate(Uint1 program_number,
   const BlastHitSavingOptions* options, Blast_Message* *blast_msg)
{
	if (options == NULL)
		return 1;

	if (options->hitlist_size < 1 || options->prelim_hitlist_size < 1)
	{
		Int4 code=1;
		Int4 subcode=1;
		Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR, code, subcode, 
                         "No hits are being saved");
		return (Int2) code;
	}

	if (options->expect_value <= 0.0 && options->cutoff_score <= 0)
	{
		Int4 code=2;
		Int4 subcode=1;
		Blast_MessageWrite(blast_msg, BLAST_SEV_ERROR, code, subcode, 
         "expect value or cutoff score must be greater than zero");
		return (Int2) code;
	}	

	return 0;
}


BlastHitSavingParameters*
BlastHitSavingParametersFree(BlastHitSavingParameters* parmameters)

{
  sfree(parmameters);
  return NULL;
}


Int2
BlastHitSavingParametersNew(Uint1 program_number, 
   const BlastHitSavingOptions* options, 
   const BlastExtensionParameters* ext_params, 
   BlastScoreBlk* sbp, BlastQueryInfo* query_info, 
   BlastHitSavingParameters* *parameters)
{
   Boolean gapped_calculation = TRUE;
   Int2 status = 0;
   BlastHitSavingParameters* params;

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

   ASSERT(options);
   ASSERT(sbp);
   
   if (!sbp->kbp_gap)
      gapped_calculation = FALSE;

   /* If parameters have not yet been created, allocate and fill all
      parameters that are constant throughout the search */
   *parameters = params = (BlastHitSavingParameters*) 
      calloc(1, sizeof(BlastHitSavingParameters));

   if (params == NULL)
      return 1;

   params->options = (BlastHitSavingOptions *) options;

   /* If sum statistics use is forced by the options, 
      set it in the paramters */
   params->do_sum_stats = options->do_sum_stats;
   /* Sum statistics is used anyway for all ungapped searches and all 
      translated gapped searches (except RPS translated searches) */
   if (!gapped_calculation || 
       (program_number != blast_type_blastn && 
        program_number != blast_type_blastp &&
        program_number != blast_type_rpsblast &&
        program_number != blast_type_rpstblastn))
      params->do_sum_stats = TRUE;
   if (program_number == blast_type_blastn || !gapped_calculation) {
      params->gap_prob = BLAST_GAP_PROB;
      params->gap_decay_rate = BLAST_GAP_DECAY_RATE;
   } else {
      params->gap_prob = BLAST_GAP_PROB_GAPPED;
      params->gap_decay_rate = BLAST_GAP_DECAY_RATE_GAPPED;
   }
   params->gap_size = BLAST_GAP_SIZE;
   params->cutoff_big_gap = 0;

   status = BlastHitSavingParametersUpdate(program_number, 
               ext_params, sbp, query_info, params);
   return status;
}

Int2
BlastHitSavingParametersUpdate(Uint1 program_number, 
   const BlastExtensionParameters* ext_params, 
   BlastScoreBlk* sbp, BlastQueryInfo* query_info, 
   BlastHitSavingParameters* params)
{
   BlastHitSavingOptions* options;
   Blast_KarlinBlk* kbp;
   double evalue;
   double scale_factor = sbp->scale_factor;

   ASSERT(params);
   ASSERT(query_info);
   ASSERT(ext_params);

   options = params->options;
   evalue = options->expect_value;

   /* Scoring options are not available here, but we can determine whether
      this is a gapped or ungapped search by checking whether gapped
      Karlin blocks have been set. */
   if (sbp->kbp_gap) {
      kbp = sbp->kbp_gap[query_info->first_context];
   } else {
      kbp = sbp->kbp[query_info->first_context];
   }

   /* Calculate cutoffs based on the current effective lengths information */
   if (options->cutoff_score > 0) {
      params->cutoff_score = options->cutoff_score * (Int4) sbp->scale_factor;
   } else if (!options->phi_align) {
      Int4 context = query_info->first_context;
      Int8 searchsp = query_info->eff_searchsp_array[context];

      /* translated RPS searches must scale the search space down */
      if (program_number == blast_type_rpstblastn)
         searchsp = searchsp / NUM_FRAMES;

      params->cutoff_score = 0;
      BLAST_Cutoffs(&params->cutoff_score, &evalue, kbp, searchsp, FALSE, 0);
      params->cutoff_score *= (Int4) scale_factor;

      /* When sum statistics is used, all HSPs above the gap trigger 
         cutoff are saved until the sum statistics is applied to potentially
         link them with other HSPs and improve their e-values. 
         However this does not apply to the ungapped search! */
      if (params->do_sum_stats) {
         params->cutoff_score = 
            MIN(params->cutoff_score, ext_params->gap_trigger);
      }
   } else {
      params->cutoff_score = 0;
   }
   
   params->cutoff_small_gap = 
      MIN(params->cutoff_score, ext_params->gap_trigger);
      
   return 0;
}

Int2 PSIBlastOptionsNew(PSIBlastOptions** psi_options)
{
   PSIBlastOptions* options;
   if (!psi_options)
      return 0;
   options = (PSIBlastOptions*)calloc(1, sizeof(PSIBlastOptions));
   *psi_options = options;
   options->inclusion_ethresh = PSI_INCLUSION_ETHRESH;
   options->pseudo_count = PSI_PSEUDO_COUNT_CONST;
   
   return 0;
}

PSIBlastOptions* PSIBlastOptionsFree(PSIBlastOptions* psi_options)
{
   sfree(psi_options);
   return NULL;
}

Int2 BlastDatabaseOptionsNew(BlastDatabaseOptions** db_options)
{
   BlastDatabaseOptions* options = (BlastDatabaseOptions*)
      calloc(1, sizeof(BlastDatabaseOptions));

   options->genetic_code = BLAST_GENETIC_CODE;
   *db_options = options;

   return 0;
}

BlastDatabaseOptions* 
BlastDatabaseOptionsFree(BlastDatabaseOptions* db_options)
{
   sfree(db_options->gen_code_string);
   sfree(db_options);
   return NULL;
}

Int2 BLAST_InitDefaultOptions(Uint1 program_number,
   LookupTableOptions** lookup_options,
   QuerySetUpOptions** query_setup_options, 
   BlastInitialWordOptions** word_options,
   BlastExtensionOptions** ext_options,
   BlastHitSavingOptions** hit_options,
   BlastScoringOptions** score_options,
   BlastEffectiveLengthsOptions** eff_len_options,
   PSIBlastOptions** psi_options,
   BlastDatabaseOptions** db_options)
{
   Int2 status;

   if ((status = LookupTableOptionsNew(program_number, lookup_options)))
      return status;

   if ((status=BlastQuerySetUpOptionsNew(query_setup_options)))
      return status;

   if ((status=BlastInitialWordOptionsNew(program_number, word_options)))
      return status;

   if ((status = BlastExtensionOptionsNew(program_number, ext_options)))
      return status;

   if ((status=BlastHitSavingOptionsNew(program_number, hit_options)))
      return status;

   if ((status=BlastScoringOptionsNew(program_number, score_options)))
      return status;

   if ((status=BlastEffectiveLengthsOptionsNew(eff_len_options)))
      return status;
   
   if ((status=PSIBlastOptionsNew(psi_options)))
      return status;

   if ((status=BlastDatabaseOptionsNew(db_options)))
      return status;

   return 0;

}

Int2 BLAST_ValidateOptions(Uint1 program_number,
                           const BlastExtensionOptions* ext_options,
                           const BlastScoringOptions* score_options, 
                           const LookupTableOptions* lookup_options, 
                           const BlastHitSavingOptions* hit_options,
                           Blast_Message* *blast_msg)
{
   Int2 status = 0;

   if ((status = BlastExtensionOptionsValidate(program_number, ext_options,
                                               blast_msg)) != 0)
       return status;
   if ((status = BlastScoringOptionsValidate(program_number, score_options,
                                               blast_msg)) != 0)
       return status;
   if ((status = LookupTableOptionsValidate(program_number, 
                    lookup_options, blast_msg)) != 0)   
       return status;
   if ((status = BlastHitSavingOptionsValidate(program_number, hit_options,
                                               blast_msg)) != 0)
       return status;
   return status;
}

/** machine epsilon assumed by CalculateLinkHSPCutoffs */
#define MY_EPS 1.0e-9

void
CalculateLinkHSPCutoffs(Uint1 program, BlastQueryInfo* query_info, 
   BlastScoreBlk* sbp, BlastHitSavingParameters* hit_params, 
   BlastExtensionParameters* ext_params,
   Int8 db_length, Int4 subject_length)
{
	double gap_prob, gap_decay_rate, x_variable, y_variable;
	Blast_KarlinBlk* kbp;
	Int4 expected_length, gap_size, query_length;
	Int8 search_sp;
        Boolean translated_subject = (program == blast_type_tblastn || 
                                 program == blast_type_rpstblastn || 
                                 program == blast_type_tblastx);

	/* Do this for the first context, should this be changed?? */
	kbp = sbp->kbp[query_info->first_context];
	gap_size = hit_params->gap_size;
	gap_prob = hit_params->gap_prob;
	gap_decay_rate = hit_params->gap_decay_rate;
   /* Use average query length */
	query_length = query_info->context_offsets[query_info->last_context+1] /
      (query_info->last_context + 1);

   if (translated_subject) {
      /* Lengths in subsequent calculations should be on the protein scale */
      subject_length /= CODON_LENGTH;
      db_length /= CODON_LENGTH;
   }

	/* Subtract off the expected score. */
   expected_length = BLAST_Nint(log(kbp->K*((double) query_length)*
                                    ((double) subject_length))/(kbp->H));
   query_length = query_length - expected_length;

   subject_length = subject_length - expected_length;
   query_length = MAX(query_length, 1);
   subject_length = MAX(subject_length, 1);

   /* If this is a database search, use database length, else the single 
      subject sequence length */
   if (db_length > subject_length) {
      y_variable = log((double) (db_length)/(double) subject_length)*(kbp->K)/
         (gap_decay_rate);
   } else {
      y_variable = log((double) (subject_length + expected_length)/
                       (double) subject_length)*(kbp->K)/(gap_decay_rate);
   }
   search_sp = ((Int8) query_length)* ((Int8) subject_length);
   x_variable = 0.25*y_variable*((double) search_sp);

   /* To use "small" gaps the query and subject must be "large" compared to
      the gap size. If small gaps may be used, then the cutoff values must be
      adjusted for the "bayesian" possibility that both large and small gaps 
      are being checked for. */

   if (search_sp > 8*gap_size*gap_size) {
      x_variable /= (1.0 - gap_prob + MY_EPS);
      hit_params->cutoff_big_gap = 
         (Int4) floor((log(x_variable)/kbp->Lambda)) + 1;
      x_variable = y_variable*(gap_size*gap_size);
      x_variable /= (gap_prob + MY_EPS);
      hit_params->cutoff_small_gap = 
         MAX(ext_params->gap_trigger, 
             (Int4) floor((log(x_variable)/kbp->Lambda)) + 1);

      hit_params->ignore_small_gaps = FALSE;
   } else {
      hit_params->cutoff_big_gap = 
         (Int4) floor((log(x_variable)/kbp->Lambda)) + 1;
      hit_params->cutoff_small_gap = hit_params->cutoff_big_gap;
      hit_params->ignore_small_gaps = TRUE;
   }	

   hit_params->cutoff_big_gap *= (Int4)sbp->scale_factor;