blast_traceback.c

ncbi源码

            /* Score is below threshold */
            gap_align->edit_block = GapEditBlockDelete(gap_align->edit_block);
            hsp_array[index] = Blast_HSPFree(hsp);
         }
      } else { 
         /* Contained within another HSP, delete. */
         hsp_array[index] = Blast_HSPFree(hsp);
      }
   } /* End loop on HSPs */

    if (translation_buffer) {
        sfree(translation_buffer);
    }
    if (frame_offsets) {
        sfree(frame_offsets);
    }
    /* Now try to detect simular alignments */

    BLASTCheckHSPInclusion(hsp_array, hsp_list->hspcnt, k_is_ooframe);
    Blast_HSPListPurgeNullHSPs(hsp_list);
    
    /* Relink and rereap the HSP list, if needed. */

    if (program_number == blast_type_blastx ||
        program_number == blast_type_tblastn ||
        program_number == blast_type_rpstblastn) {
        
        if (hit_params->do_sum_stats == TRUE) {
           BLAST_LinkHsps(program_number, hsp_list, query_info, subject_blk,
                         sbp, hit_params, score_options->gapped_calculation);
        } else if (hit_options->phi_align) {
           Blast_HSPListPHIGetEvalues(hsp_list, sbp);
        } else {
           Blast_HSPListGetEvalues(program_number, query_info, hsp_list, 
                                      score_options->gapped_calculation, sbp);
        }
        
        Blast_HSPListReapByEvalue(hsp_list, hit_options);
    }
    
    qsort(hsp_array, hsp_list->hspcnt, sizeof(BlastHSP*), score_compare_hsps);
    
    /* Remove extra HSPs if there is a user proveded limit on the number 
       of HSPs per database sequence */
    if (hit_options->hsp_num_max > 0 && 
        hit_options->hsp_num_max < hsp_list->hspcnt) {
       for (index=hit_options->hsp_num_max; index<hsp_list->hspcnt; ++index) {
          hsp_array[index] = Blast_HSPFree(hsp_array[index]);
       }
       hsp_list->hspcnt = hit_options->hsp_num_max;
    }

    return 0;
}

static Uint1 Blast_TracebackGetEncoding(Uint1 program_number) 
{
   Uint1 encoding;

   switch (program_number) {
   case blast_type_blastn:
      encoding = BLASTNA_ENCODING;
      break;
   case blast_type_blastp:
   case blast_type_rpsblast:
      encoding = BLASTP_ENCODING;
      break;
   case blast_type_blastx:
      encoding = BLASTP_ENCODING;
      break;
   case blast_type_tblastn:
   case blast_type_rpstblastn:
      encoding = NCBI4NA_ENCODING;
      break;
   case blast_type_tblastx:
      encoding = NCBI4NA_ENCODING;
      break;
   default:
      encoding = ERROR_ENCODING;
      break;
   }
   return encoding;
}

/** Delete extra subject sequences hits, if after-traceback hit list size is
 * smaller than preliminary hit list size.
 * @param results All results after traceback, assumed already sorted by best
 *                e-value [in] [out]
 * @param hitlist_size Final hit list size [in]
 */
static void 
BlastPruneExtraHits(BlastHSPResults* results, Int4 hitlist_size)
{
   Int4 query_index, subject_index;
   BlastHitList* hit_list;

   for (query_index = 0; query_index < results->num_queries; ++query_index) {
      if (!(hit_list = results->hitlist_array[query_index]))
         continue;
      for (subject_index = hitlist_size;
           subject_index < hit_list->hsplist_count; ++subject_index) {
         hit_list->hsplist_array[subject_index] = 
            Blast_HSPListFree(hit_list->hsplist_array[subject_index]);
      }
      hit_list->hsplist_count = MIN(hit_list->hsplist_count, hitlist_size);
   }
}

Int2 BLAST_ComputeTraceback(Uint1 program_number, BlastHSPResults* results, 
        BLAST_SequenceBlk* query, BlastQueryInfo* query_info, 
        const BlastSeqSrc* seq_src, BlastGapAlignStruct* gap_align,
        BlastScoringParameters* score_params,
        const BlastExtensionParameters* ext_params,
        BlastHitSavingParameters* hit_params,
        BlastEffectiveLengthsParameters* eff_len_params,
        const BlastDatabaseOptions* db_options,
        const PSIBlastOptions* psi_options)
{
   Int2 status = 0;
   Int4 query_index, subject_index;
   BlastHitList* hit_list;
   BlastHSPList* hsp_list;
   BlastScoreBlk* sbp;
   Uint1 encoding;
   GetSeqArg seq_arg;
   
   if (!results || !query_info || !seq_src) {
      return 0;
   }
   
   /* Set the raw X-dropoff value for the final gapped extension with 
      traceback */
   gap_align->gap_x_dropoff = ext_params->gap_x_dropoff_final;

   sbp = gap_align->sbp;
   
   encoding = Blast_TracebackGetEncoding(program_number);
   memset((void*) &seq_arg, 0, sizeof(seq_arg));

   for (query_index = 0; query_index < results->num_queries; ++query_index) {
      hit_list = results->hitlist_array[query_index];

      if (!hit_list)
         continue;

     if (program_number == blast_type_blastp && 
         (ext_params->options->compositionBasedStats == TRUE || 
            ext_params->options->eTbackExt == eSmithWatermanTbck))
     {
         Kappa_RedoAlignmentCore(query, query_info, sbp, hit_list, seq_src, 
           score_params, ext_params, hit_params, psi_options); 
     }
     else
     {
      for (subject_index = 0; subject_index < hit_list->hsplist_count;
           ++subject_index) {
         hsp_list = hit_list->hsplist_array[subject_index];
         if (!hsp_list)
            continue;

         if (!hsp_list->traceback_done) {
            seq_arg.oid = hsp_list->oid;
            seq_arg.encoding = encoding;
            BlastSequenceBlkClean(seq_arg.seq);
            if (BLASTSeqSrcGetSequence(seq_src, (void*) &seq_arg) < 0)
                continue;

            if (BLASTSeqSrcGetTotLen(seq_src) == 0) {
               /* This is not a database search, so effective search spaces
                  need to be recalculated based on this subject sequence 
                  length */
               if ((status = BLAST_OneSubjectUpdateParameters(program_number, 
                                seq_arg.seq->length, score_params->options, 
                                query_info, sbp, ext_params, hit_params, 
                                NULL, eff_len_params)) != 0)
                  return status;
            }

            Blast_TracebackFromHSPList(program_number, hsp_list, query, 
               seq_arg.seq, query_info, gap_align, sbp, score_params, 
               ext_params->options, hit_params, db_options->gen_code_string);
            BLASTSeqSrcRetSequence(seq_src, (void*)&seq_arg);
         }
      }
     }
   }

   /* Re-sort the hit lists according to their best e-values, because
      they could have changed. Only do this for a database search. */
   if (BLASTSeqSrcGetTotLen(seq_src) > 0)
      Blast_HSPResultsSortByEvalue(results);

   /* Eliminate extra hits from results, if preliminary hit list size is larger
      than the final hit list size */
   if (hit_params->options->hitlist_size < 
       hit_params->options->prelim_hitlist_size)
      BlastPruneExtraHits(results, hit_params->options->hitlist_size);

   BlastSequenceBlkFree(seq_arg.seq);

   return status;
}

#define SWAP(a, b) {tmp = (a); (a) = (b); (b) = tmp; }

static void 
RPSUpdateTraceback(BlastHSP *hsp)
{
   Int4 tmp;
   GapEditBlock *gap_info = hsp->gap_info;
   GapEditScript *esp;

   if (gap_info == NULL)
      return;

   SWAP(gap_info->start1, gap_info->start2);
   SWAP(gap_info->length1, gap_info->length2);
   SWAP(gap_info->original_length1, gap_info->original_length2);
   SWAP(gap_info->frame1, gap_info->frame2);
   SWAP(gap_info->translate1, gap_info->translate2);

   esp = gap_info->esp;
   while (esp != NULL) {
      if (esp->op_type == eGapAlignIns)
          esp->op_type = eGapAlignDel;
      else if (esp->op_type == eGapAlignDel)
          esp->op_type = eGapAlignIns;

      esp = esp->next;
   }
}

static void 
RPSUpdateHSPList(BlastHSPList *hsplist)
{
   Int4 i;
   BlastHSP **hsp;
   BlastSeg tmp;

   hsp = hsplist->hsp_array;
   for (i = 0; i < hsplist->hspcnt; i++) {

      /* switch query and subject offsets (which are
         already in local coordinates) */

      tmp = hsp[i]->query;
      hsp[i]->query = hsp[i]->subject;
      hsp[i]->subject = tmp;

      /* Change the traceback information to reflect the
         query and subject sequences getting switched */

      RPSUpdateTraceback(hsp[i]);
   }
}

#define RPS_K_MULT 1.2

Int2 BLAST_RPSTraceback(Uint1 program_number, 
        BlastHSPResults* results, 
        BLAST_SequenceBlk* concat_db, BlastQueryInfo* concat_db_info, 
        BLAST_SequenceBlk* query, BlastQueryInfo* query_info, 
        BlastGapAlignStruct* gap_align,
        const BlastScoringParameters* score_params,
        const BlastExtensionParameters* ext_params,
        BlastHitSavingParameters* hit_params,
        const BlastDatabaseOptions* db_options,
        const double* karlin_k)
{
   Int2 status = 0;
   Int4 i;
   BlastHitList* hit_list;
   BlastHSPList* hsp_list;
   BlastScoreBlk* sbp;
   Int4 **orig_pssm;
   
   if (!results || !concat_db_info || !concat_db) {
      return 0;
   }
   
   /* Set the raw X-dropoff value for the final gapped extension with 
      traceback */
   gap_align->gap_x_dropoff = ext_params->gap_x_dropoff_final;

   sbp = gap_align->sbp;
   orig_pssm = gap_align->sbp->posMatrix;

   hit_list = results->hitlist_array[0];
   if (!hit_list)
      return 0;

   /* for translated searches, the traceback code calculates
      E values *after* the scaling factor has been removed from
      the alignment scores. Thus, lambda must not be pre-scaled
      for a translated search */

   if (program_number != blast_type_rpstblastn)
      sbp->kbp_gap[0]->Lambda /= score_params->scale_factor;

   for (i = 0; i < hit_list->hsplist_count; i++) {
      hsp_list = hit_list->hsplist_array[i];
      if (!hsp_list)
         continue;

      if (!hsp_list->traceback_done) {

         Int4 offsets[2];
         BLAST_SequenceBlk one_db_seq;
         BlastQueryInfo one_db_seq_info;
         Int4 *db_seq_start;

         /* pick out one of the sequences from the concatenated
            DB (given by the OID of this HSPList). The sequence
            size does not include the trailing NULL */

         db_seq_start = &concat_db_info->context_offsets[hsp_list->oid];
         memset(&one_db_seq, 0, sizeof(one_db_seq));
         one_db_seq.sequence = NULL;
         one_db_seq.length = db_seq_start[1] - db_seq_start[0] - 1;

         /* Set up the QueryInfo structure for this sequence. The
            trailing NULL must be added back */

         offsets[0] = 0;
         offsets[1] = one_db_seq.length + 1;

         memset(&one_db_seq_info, 0, sizeof(one_db_seq_info));
         one_db_seq_info.first_context = 0;
         one_db_seq_info.last_context = 0;
         one_db_seq_info.num_queries = 1;
         one_db_seq_info.context_offsets = &offsets[0];
         one_db_seq_info.eff_searchsp_array = query_info->eff_searchsp_array;

         /* Update the statistics for this database sequence
            (if not a translated search) */

         if (program_number == blast_type_rpstblastn) {
            sbp->posMatrix = orig_pssm + db_seq_start[0];
         }
         else {
            /* replace the PSSM and the Karlin values
               for this DB sequence. */

            sbp->posMatrix = RPSCalculatePSSM(score_params->scale_factor,
                        query->length, query->sequence, one_db_seq.length,
                        orig_pssm + db_seq_start[0]);
            if (sbp->posMatrix == NULL)
               return -1;

            sbp->kbp_gap[0]->K = RPS_K_MULT * karlin_k[hsp_list->oid];
            sbp->kbp_gap[0]->logK = log(RPS_K_MULT * karlin_k[hsp_list->oid]);
         }

         /* compute the traceback information and calculate E values
            for all HSPs in the list */

         Blast_TracebackFromHSPList(program_number, hsp_list, &one_db_seq, 
            query, &one_db_seq_info, gap_align, sbp, score_params, 
            ext_params->options, hit_params, db_options->gen_code_string);

         if (program_number != blast_type_rpstblastn)
            _PSIDeallocateMatrix((void**)sbp->posMatrix, one_db_seq.length+1);
      }

      /* Revert query and subject to their traditional meanings. 
         This involves switching the offsets around and reversing
         any traceback information */

      RPSUpdateHSPList(hsp_list);
   }

   /* restore input data */
   if (program_number != blast_type_rpstblastn)
      sbp->kbp_gap[0]->Lambda *= score_params->scale_factor;

   gap_align->sbp->posMatrix = orig_pssm;
   return status;
}