blast_hits.c

ncbi源码

         new_segment1 = new_segment1->next;
         num1++;
      } else {
         new_segment2 = new_segment2->next;
         num2++;
      }
   }
   
   sfree(segments1);
   sfree(segments2);

   if (intersection_found) {
      esp = NULL;
      for (index = 0; index < num1-1; index++)
         esp1 = esp1->next;
      for (index = 0; index < num2-1; index++) {
         esp = esp2;
         esp2 = esp2->next;
      }
      if (intersection_found < 3) {
         if (num1 > 0)
            esp1 = esp1->next;
         if (num2 > 0) {
            esp = esp2;
            esp2 = esp2->next;
         }
      }
      switch (intersection_found) {
      case 1:
         esp1->num = dist;
         esp1->next = esp2->next;
         esp2->next = NULL;
         break;
      case 2:
         esp1->num = dist;
         esp2->num = next_dist;
         esp1->next = esp2;
         if (esp)
            esp->next = NULL;
         break;
      case 3:
         esp1->num += dist;
         esp2->op_type = op_type;
         esp2->num = next_dist;
         esp1->next = esp2;
         if (esp)
            esp->next = NULL;
         break;
      default: break;
      }
      hsp1->query.end = hsp2->query.end;
      hsp1->subject.end = hsp2->subject.end;
      hsp1->query.length = hsp1->query.end - hsp1->query.offset;
      hsp1->subject.length = hsp1->subject.end - hsp1->subject.offset;
   }

   return (Boolean) intersection_found;
}

/********************************************************************************
          Functions manipulating BlastHSPList's
********************************************************************************/

BlastHSPList* Blast_HSPListFree(BlastHSPList* hsp_list)
{
   Int4 index;

   if (!hsp_list)
      return hsp_list;

   for (index = 0; index < hsp_list->hspcnt; ++index) {
      Blast_HSPFree(hsp_list->hsp_array[index]);
   }
   sfree(hsp_list->hsp_array);

   sfree(hsp_list);
   return NULL;
}

BlastHSPList* Blast_HSPListNew(Int4 hsp_max)
{
   BlastHSPList* hsp_list = (BlastHSPList*) calloc(1, sizeof(BlastHSPList));
   const Int4 k_default_allocated=100;

   /* hsp_max is max number of HSP's ever allowed, INT4_MAX taken as infinity. */
   hsp_list->hsp_max = INT4_MAX; 
   if (hsp_max > 0)
      hsp_list->hsp_max = hsp_max;

   hsp_list->allocated = MIN(k_default_allocated, hsp_list->hsp_max);

   hsp_list->hsp_array = (BlastHSP**) 
      calloc(hsp_list->allocated, sizeof(BlastHSP*));

   return hsp_list;
}

/** Duplicate HSPList's contents, copying only pointers to individual HSPs */
static BlastHSPList* Blast_HSPListDup(BlastHSPList* hsp_list)
{
   BlastHSPList* new_hsp_list = (BlastHSPList*) 
      BlastMemDup(hsp_list, sizeof(BlastHSPList));
   new_hsp_list->hsp_array = (BlastHSP**) 
      BlastMemDup(hsp_list->hsp_array, hsp_list->allocated*sizeof(BlastHSP*));
   new_hsp_list->allocated = hsp_list->allocated;

   return new_hsp_list;
}

/* Comments in blast_hits.h
 */
Int2 
Blast_HSPListSaveHSP(BlastHSPList* hsp_list, BlastHSP* new_hsp)
{
   BlastHSP** hsp_array;
   Int4 highscore, lowscore, score = 0;
   Int4 hspcnt, new_index, old_index;
   Int4 hsp_allocated; /* how many hsps are in the array. */

   hspcnt = hsp_list->hspcnt;
   hsp_allocated = hsp_list->allocated;
   
   /* Check if list is already full, then reallocate. */
   if (hspcnt >= hsp_allocated-1 && hsp_list->do_not_reallocate == FALSE)
   {
      Int4 new_allocated = MIN(2*hsp_list->allocated, hsp_list->hsp_max);
      if (new_allocated <= hsp_list->hsp_max) {
         hsp_array = (BlastHSP**)
            realloc(hsp_list->hsp_array, new_allocated*sizeof(BlastHSP*));
         if (hsp_array == NULL)
         {
#ifdef ERR_POST_EX_DEFINED
            ErrPostEx(SEV_WARNING, 0, 0, 
               "UNABLE to reallocate in Blast_SaveHsp,"
               " continuing with fixed array of %ld HSP's", 
                      (long) hsp_allocated);
#endif
            hsp_list->do_not_reallocate = TRUE; 
         } else {
            hsp_list->hsp_array = hsp_array;
            hsp_list->allocated = new_allocated;
            hsp_allocated = new_allocated;
         }
      } else {
         hsp_list->do_not_reallocate = TRUE; 
      }
   }

   hsp_array = hsp_list->hsp_array;

   /* If we are saving ALL HSP's, simply save and sort later. */
   if (hsp_list->do_not_reallocate == FALSE)
   {
      hsp_array[hsp_list->hspcnt] = new_hsp;
      (hsp_list->hspcnt)++;
      return 0;
   }

   /* Use a binary search to insert the HSP. */
   score = new_hsp->score;
    

   if (hspcnt != 0) {
      highscore = hsp_array[0]->score;
      lowscore = hsp_array[hspcnt-1]->score;
   } else {
      highscore = 0;
      lowscore = 0;
   }

   if (score >= highscore) {
      new_index = 0;
   } else if (score <= lowscore) {
      new_index = hspcnt;
   } else {
      const Int4 k_MaxIter=30;
      Int4 index;
      Int4 low_index = 0;
      Int4 high_index = hspcnt-1;

      new_index = (low_index+high_index)/2;
      old_index = new_index;
      
      for (index=0; index<k_MaxIter; index++)
      {
         if (score > hsp_array[new_index]->score)
            high_index = new_index;
         else
            low_index = new_index;
         new_index = (low_index+high_index)/2;
         if (new_index == old_index) { 
            /* Perform this check as new_index get rounded DOWN above.*/
            if (score < hsp_array[new_index]->score)
               new_index++;
            break;
         }
         old_index = new_index;
      }
   }

   if (hspcnt >= hsp_allocated-1) {
      if (new_index >= hspcnt) { 
         /* this HSP is less significant than others on a full list.*/
         sfree(new_hsp);
         return 0;
      } else { /* Delete the last HSP on the list. */
         hspcnt = hsp_list->hspcnt--;
         sfree(hsp_array[hspcnt-1]);
      }
   }
   hsp_list->hspcnt++;
   memmove((hsp_array+new_index+1), (hsp_array+new_index), (hspcnt-new_index)*sizeof(hsp_array[0]));
   hsp_array[new_index] = new_hsp;
   
   return 0;
}

void 
Blast_HSPListSetFrames(Uint1 program_number, BlastHSPList* hsp_list, 
                 Boolean is_ooframe)
{
   BlastHSP* hsp;
   Int4 index;

   if (!hsp_list)
      return;

   for (index=0; index<hsp_list->hspcnt; index++) {
      hsp = hsp_list->hsp_array[index];
      if (is_ooframe && program_number == blast_type_blastx) {
         /* Query offset is in mixed-frame coordinates */
         hsp->query.frame = hsp->query.offset % CODON_LENGTH + 1;
         if ((hsp->context % NUM_FRAMES) >= CODON_LENGTH)
            hsp->query.frame = -hsp->query.frame;
      } else {
         hsp->query.frame = BLAST_ContextToFrame(program_number, hsp->context);
      }
      
      /* Correct offsets in the edit block too */
      if (hsp->gap_info) {
         hsp->gap_info->frame1 = hsp->query.frame;
         hsp->gap_info->frame2 = hsp->subject.frame;
      }
   }
}

Int2 Blast_HSPListGetEvalues(Uint1 program, BlastQueryInfo* query_info,
        BlastHSPList* hsp_list, Boolean gapped_calculation, 
        BlastScoreBlk* sbp)
{
   BlastHSP* hsp;
   BlastHSP** hsp_array;
   Blast_KarlinBlk** kbp;
   Int4 hsp_cnt;
   Int4 index;
   
   if (hsp_list == NULL)
      return 1;
   
   if (gapped_calculation)
      kbp = sbp->kbp_gap_std;
   else
      kbp = sbp->kbp_std;
   
   hsp_cnt = hsp_list->hspcnt;
   hsp_array = hsp_list->hsp_array;
   for (index=0; index<hsp_cnt; index++) {
      hsp = hsp_array[index];

      ASSERT(hsp != NULL);

      /* Get effective search space from the score block, or from the 
         query information block, in order of preference */
      if (sbp->effective_search_sp) {
         hsp->evalue = BLAST_KarlinStoE_simple(hsp->score, kbp[hsp->context],
                          sbp->effective_search_sp);
      } else {
         hsp->evalue = 
            BLAST_KarlinStoE_simple(hsp->score, kbp[hsp->context],
               query_info->eff_searchsp_array[hsp->context]);
      }
   }
   
   return 0;
}

void Blast_HSPListPHIGetEvalues(BlastHSPList* hsp_list, BlastScoreBlk* sbp)
{
   Int4 index;
   BlastHSP* hsp;

   for (index = 0; index < hsp_list->hspcnt; ++index) {
      hsp = hsp_list->hsp_array[index];
      Blast_HSPPHIGetEvalue(hsp, sbp);
   }
}

Int2 Blast_HSPListReapByEvalue(BlastHSPList* hsp_list, 
        BlastHitSavingOptions* hit_options)
{
   BlastHSP* hsp;
   BlastHSP** hsp_array;
   Int4 hsp_cnt = 0;
   Int4 index;
   double cutoff;
   
   if (hsp_list == NULL)
      return 1;

   cutoff = hit_options->expect_value;

   hsp_array = hsp_list->hsp_array;
   for (index = 0; index < hsp_list->hspcnt; index++) {
      hsp = hsp_array[index];

      ASSERT(hsp != NULL);
      
      if (hsp->evalue > cutoff) {
         hsp_array[index] = Blast_HSPFree(hsp_array[index]);
      } else {
         if (index > hsp_cnt)
            hsp_array[hsp_cnt] = hsp_array[index];
         hsp_cnt++;
      }
   }
      
   hsp_list->hspcnt = hsp_cnt;

   return 0;
}

/* See description in blast_hits.h */

Int2
Blast_HSPListPurgeNullHSPs(BlastHSPList* hsp_list)

{
        Int4 index, index1; /* loop indices. */
        Int4 hspcnt; /* total number of HSP's to iterate over. */
        BlastHSP** hsp_array;  /* hsp_array to purge. */

	if (hsp_list == NULL || hsp_list->hspcnt == 0)
		return 0;

        hsp_array = hsp_list->hsp_array;
        hspcnt = hsp_list->hspcnt;

	index1 = 0;
	for (index=0; index<hspcnt; index++)
	{
		if (hsp_array[index] != NULL)
		{
			hsp_array[index1] = hsp_array[index];
			index1++;
		}
	}

	for (index=index1; index<hspcnt; index++)
	{
		hsp_array[index] = NULL;
	}

	hsp_list->hspcnt = index1;

        return 0;
}

/** Are the two HSPs within a given diagonal distance of each other? */
#define MB_HSP_CLOSE(q1, q2, s1, s2, c) (ABS(((q1)-(s1)) - ((q2)-(s2))) < (c))
#define MIN_DIAG_DIST 60

/** Sort the HSPs in an HSP list by diagonal and remove redundant HSPs */
static Int2
Blast_HSPListUniqSort(BlastHSPList* hsp_list)
{
   Int4 index, new_hspcnt, index1, q_off, s_off, q_end, s_end, index2;
   BlastHSP** hsp_array = hsp_list->hsp_array;
   Boolean shift_needed = FALSE;
   Int4 context;
   double evalue;

   qsort(hsp_array, hsp_list->hspcnt, sizeof(BlastHSP*), 
         diag_uniq_compare_hsps);
   /* Delete all HSPs that were flagged in qsort */
   for (index = 0; index < hsp_list->hspcnt; ++index) {
      if (hsp_array[index]->score == 0) {
         hsp_array[index] = Blast_HSPFree(hsp_array[index]);
      }
   }      
   /* Move all nulled out HSPs to the end */
   qsort(hsp_array, hsp_list->hspcnt, sizeof(BlastHSP*),
         null_compare_hsps);

   for (index=1, new_hspcnt=0; index<hsp_list->hspcnt; index++) {
      if (!hsp_array[index])
         break;

      q_off = hsp_array[index]->query.offset;
      s_off = hsp_array[index]->subject.offset;
      q_end = hsp_array[index]->query.end;
      s_end = hsp_array[index]->subject.end;
      evalue = hsp_array[index]->evalue;
      context = hsp_array[index]->context;
      for (index1 = new_hspcnt; index1 >= 0 && 
           hsp_array[index1]->context == context && new_hspcnt-index1 < 10 && 
           MB_HSP_CLOSE(q_off, hsp_array[index1]->query.offset,
                        s_off, hsp_array[index1]->subject.offset, 
                        MIN_DIAG_DIST);