link_hsps.c

ncbi源码

	} /* end for frame_index ... */

   sfree(hp_frame_start);
   sfree(hp_frame_number);
   sfree(hp_start->hsp);
   sfree(hp_start);

   if (translated_query) {
      qsort(link_hsp_array,total_number_of_hsps,sizeof(LinkHSPStruct*), 
            rev_compare_hsps_transl);
      qsort(link_hsp_array, total_number_of_hsps,sizeof(LinkHSPStruct*), 
            fwd_compare_hsps_transl);
   } else {
      qsort(link_hsp_array,total_number_of_hsps,sizeof(LinkHSPStruct*), 
            rev_compare_hsps);
      qsort(link_hsp_array, total_number_of_hsps,sizeof(LinkHSPStruct*), 
            fwd_compare_hsps);
   }

   /* Sort by starting position. */
   

	for (index=0, last_hsp=NULL;index<total_number_of_hsps; index++) 
	{
		H = link_hsp_array[index];
		H->prev = NULL;
		H->next = NULL;
	}

   /* hook up the HSP's. */
	first_hsp = NULL;
	for (index=0, last_hsp=NULL;index<total_number_of_hsps; index++) 
   {
		H = link_hsp_array[index];

      /* If this is not a single piece or the start of a chain, then Skip it. */
      if (H->linked_set == TRUE && H->start_of_chain == FALSE)
			continue;

      /* If the HSP has no "link" connect the "next", otherwise follow the "link"
         chain down, connecting them with "next" and "prev". */
		if (last_hsp == NULL)
			first_hsp = H;
		H->prev = last_hsp;
		ordering_method = H->ordering_method;
		if (H->hsp_link.link[ordering_method] == NULL)
		{
         /* Grab the next HSP that is not part of a chain or the start of a chain */
         /* The "next" pointers are not hooked up yet in HSP's further down array. */
         index1=index;
         H2 = index1<(total_number_of_hsps-1) ? link_hsp_array[index1+1] : NULL;
         while (H2 && H2->linked_set == TRUE && 
                H2->start_of_chain == FALSE)
         {
            index1++;
		     	H2 = index1<(total_number_of_hsps-1) ? link_hsp_array[index1+1] : NULL;
         }
         H->next= H2;
		}
		else
		{
			/* The first one has the number of links correct. */
			num_links = H->hsp_link.num[ordering_method];
			link = H->hsp_link.link[ordering_method];
			while (link)
			{
				H->hsp->num = num_links;
				H->sumscore = H->hsp_link.sum[ordering_method];
				H->next = (LinkHSPStruct*) link;
				H->prev = last_hsp;
				last_hsp = H;
				H = H->next;
				if (H != NULL)
				    link = H->hsp_link.link[ordering_method];
				else
				    break;
			}
			/* Set these for last link in chain. */
			H->hsp->num = num_links;
			H->sumscore = H->hsp_link.sum[ordering_method];
         /* Grab the next HSP that is not part of a chain or the start of a chain */
         index1=index;
         H2 = index1<(total_number_of_hsps-1) ? link_hsp_array[index1+1] : NULL;
         while (H2 && H2->linked_set == TRUE && 
                H2->start_of_chain == FALSE)
		   {
            index1++;
            H2 = index1<(total_number_of_hsps-1) ? link_hsp_array[index1+1] : NULL;
			}
         H->next= H2;
			H->prev = last_hsp;
		}
		last_hsp = H;
	}
	
   /* The HSP's may be in a different order than they were before, 
      but first_hsp contains the first one. */
   for (index = 0, H = first_hsp; index < hsp_list->hspcnt; index++) {
      hsp_list->hsp_array[index] = H->hsp;
      /* Free the wrapper structure */
      H2 = H->next;
      sfree(H);
      H = H2;
   }
   sfree(link_hsp_array);
   sfree(lh_helper);

   return 0;
}

static void ConnectLinkHSPStructs(LinkHSPStruct** linkhsp_array, Int4 hspcnt,
                                  Uint1* subject_seq)
{
   Int4 index, index1, i;
   LinkHSPStruct* hsp;

   qsort(linkhsp_array, hspcnt, sizeof(LinkHSPStruct*), 
         sumscore_compare_hsps);

   hsp = linkhsp_array[0];
   for (index=0; index<hspcnt; hsp = hsp->next) {
      if (hsp->linked_set) {
         index1 = hsp->hsp->num;
         for (i=1; i < index1; i++, hsp = hsp->next) {
            hsp->next->hsp->evalue = hsp->hsp->evalue; 
            hsp->next->hsp->num = hsp->hsp->num;
            hsp->next->sumscore = hsp->sumscore;
            if (FindSpliceJunction(subject_seq, hsp->hsp, hsp->next->hsp)) {
               /* Kludge: ordering_method here would indicate existence of
                  splice junctions(s) */
               hsp->hsp->splice_junction++;
               hsp->next->hsp->splice_junction++;
            } else {
               hsp->hsp->splice_junction--;
               hsp->next->hsp->splice_junction--;
            }
         }
      } 
      while (++index < hspcnt)
         if (!linkhsp_array[index]->linked_set ||
             linkhsp_array[index]->start_of_chain)
            break;
      if (index == hspcnt) {
         hsp->next = NULL;
         break;
      }
      hsp->next = linkhsp_array[index];
   }
}

static Int2
new_link_hsps(Uint1 program_number, BlastHSPList* hsp_list, 
   BlastQueryInfo* query_info, BLAST_SequenceBlk* subject,
   BlastScoreBlk* sbp, const BlastHitSavingParameters* hit_params)
{
   BlastHSP** hsp_array;
   LinkHSPStruct** score_hsp_array,** offset_hsp_array,** end_hsp_array;
   LinkHSPStruct* hsp,* head_hsp,* best_hsp,* var_hsp;
   Int4 hspcnt, index, index1, i;
   double best_evalue, evalue;
   Int4 sumscore, best_sumscore = 0;
   Boolean reverse_link;
   Int4 longest_intron = hit_params->options->longest_intron;
   LinkHSPStruct** link_hsp_array;

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

   link_hsp_array = (LinkHSPStruct**) malloc(hspcnt*sizeof(LinkHSPStruct*));
   for (index = 0; index < hspcnt; ++index) {
      link_hsp_array[index] = (LinkHSPStruct*) calloc(1, sizeof(LinkHSPStruct));
      link_hsp_array[index]->hsp = hsp_array[index];
      link_hsp_array[index]->linked_set = FALSE;
      hsp_array[index]->num = 1;
      hsp_array[index]->splice_junction = FALSE;
   }

   score_hsp_array = (LinkHSPStruct**) malloc(hspcnt*sizeof(LinkHSPStruct*));
   offset_hsp_array = (LinkHSPStruct**) malloc(hspcnt*sizeof(LinkHSPStruct*));
   end_hsp_array = (LinkHSPStruct**) malloc(hspcnt*sizeof(LinkHSPStruct*));

   memcpy(score_hsp_array, link_hsp_array, hspcnt*sizeof(LinkHSPStruct*));
   memcpy(offset_hsp_array, link_hsp_array, hspcnt*sizeof(LinkHSPStruct*));
   memcpy(end_hsp_array, link_hsp_array, hspcnt*sizeof(LinkHSPStruct*));
   qsort(offset_hsp_array, hspcnt, sizeof(LinkHSPStruct*), fwd_compare_hsps);
   qsort(end_hsp_array, hspcnt, sizeof(LinkHSPStruct*), end_compare_hsps);

   qsort(score_hsp_array, hspcnt, sizeof(LinkHSPStruct*), 
            sumscore_compare_hsps);
      
   head_hsp = NULL;
   for (index = 0; index < hspcnt && score_hsp_array[index]; ) {
      if (!head_hsp) {
         while (index<hspcnt && score_hsp_array[index] && 
                score_hsp_array[index]->linked_set)
            index++;
         if (index==hspcnt)
            break;
         head_hsp = score_hsp_array[index];
      }
      best_evalue = head_hsp->hsp->evalue;
      best_hsp = NULL;
      reverse_link = FALSE;
      
      if (head_hsp->linked_set)
         for (var_hsp = head_hsp, i=1; i<head_hsp->hsp->num; 
              var_hsp = var_hsp->next, i++);
      else
         var_hsp = head_hsp;

      index1 = hsp_binary_search(offset_hsp_array, hspcnt,
                                 var_hsp->hsp->query.end - WINDOW_SIZE, TRUE);
      while (index1 < hspcnt && offset_hsp_array[index1]->hsp->query.offset < 
             var_hsp->hsp->query.end + WINDOW_SIZE) {
         hsp = offset_hsp_array[index1++];
         /* If this is already part of a linked set, disregard it */
         if (hsp == var_hsp || hsp == head_hsp || 
             (hsp->linked_set && !hsp->start_of_chain))
            continue;
         /* Check if the subject coordinates are consistent with query */
         if ((hsp->hsp->subject.offset < 
             var_hsp->hsp->subject.end - WINDOW_SIZE) || 
             (hsp->hsp->subject.offset > 
              var_hsp->hsp->subject.end + longest_intron))
            continue;
         /* Check if the e-value for the combined two HSPs is better than for
            one of them */
         if ((evalue = SumHSPEvalue(program_number, sbp, query_info, subject, 
                                    hit_params, head_hsp, hsp, &sumscore)) < 
             MIN(best_evalue, hsp->hsp->evalue)) {
            best_hsp = hsp;
            best_evalue = evalue;
            best_sumscore = sumscore;
         }
      }
      index1 = hsp_binary_search(end_hsp_array, hspcnt,
                                 head_hsp->hsp->query.offset - WINDOW_SIZE, FALSE);
      while (index1 < hspcnt && end_hsp_array[index1]->hsp->query.end < 
             head_hsp->hsp->query.offset + WINDOW_SIZE) {
         hsp = end_hsp_array[index1++];

         /* Check if the subject coordinates are consistent with query */
         if (hsp == head_hsp || 
             hsp->hsp->subject.end > head_hsp->hsp->subject.offset + WINDOW_SIZE || 
             hsp->hsp->subject.end < head_hsp->hsp->subject.offset - longest_intron)
            continue;
         if (hsp->linked_set) {
            for (var_hsp = hsp, i=1; var_hsp->start_of_chain == FALSE; 
                 var_hsp = var_hsp->prev, i++);
            if (i<var_hsp->hsp->num || var_hsp == head_hsp)
               continue;
         } else
            var_hsp = hsp;
         /* Check if the e-value for the combined two HSPs is better than for
            one of them */
         if ((evalue = SumHSPEvalue(program_number, sbp, query_info, subject, 
                          hit_params, var_hsp, head_hsp, &sumscore)) < 
             MIN(var_hsp->hsp->evalue, best_evalue)) {
            best_hsp = hsp;
            best_evalue = evalue;
            best_sumscore = sumscore;
            reverse_link = TRUE;
         }
      }
         
      /* Link these HSPs together */
      if (best_hsp != NULL) {
         if (!reverse_link) {
            head_hsp->start_of_chain = TRUE;
            head_hsp->sumscore = best_sumscore;
            head_hsp->hsp->evalue = best_evalue;
            best_hsp->start_of_chain = FALSE;
            if (head_hsp->linked_set) 
               for (var_hsp = head_hsp, i=1; i<head_hsp->hsp->num; 
                    var_hsp = var_hsp->next, i++);
            else 
               var_hsp = head_hsp;
            var_hsp->next = best_hsp;
            best_hsp->prev = var_hsp;
            head_hsp->hsp->num += best_hsp->hsp->num;
         } else {
            best_hsp->next = head_hsp;
            head_hsp->prev = best_hsp;
            if (best_hsp->linked_set) {
               for (var_hsp = best_hsp; 
                    var_hsp->start_of_chain == FALSE; 
                    var_hsp = var_hsp->prev);
            } else
                  var_hsp = best_hsp;
            var_hsp->start_of_chain = TRUE;
            var_hsp->sumscore = best_sumscore;
            var_hsp->hsp->evalue = best_evalue;
            var_hsp->hsp->num += head_hsp->hsp->num;
            head_hsp->start_of_chain = FALSE;
         }
         
         head_hsp->linked_set = best_hsp->linked_set = TRUE;
         if (reverse_link)
            head_hsp = var_hsp;
      } else {
         head_hsp = NULL;
         index++;
      }
   }
  
   ConnectLinkHSPStructs(score_hsp_array, hspcnt, 
                         subject->sequence_start);
   head_hsp = score_hsp_array[0];

   sfree(score_hsp_array);
   sfree(offset_hsp_array);
   sfree(end_hsp_array);

   /* Place HSPs in the original HSP array in their new order */
   for (index=0; head_hsp && index<hspcnt; index++) {
      hsp_list->hsp_array[index] = head_hsp->hsp;
      var_hsp = head_hsp->next;
      sfree(head_hsp);
      head_hsp = var_hsp;
   }

   sfree(link_hsp_array);
   return 0;
}

Int2 
BLAST_LinkHsps(Uint1 program_number, BlastHSPList* hsp_list, 
   BlastQueryInfo* query_info, BLAST_SequenceBlk* subject, 
   BlastScoreBlk* sbp, const BlastHitSavingParameters* hit_params,
   Boolean gapped_calculation)
{
	if (hsp_list && hsp_list->hspcnt > 0)
	{
      /* Link up the HSP's for this hsp_list. */
      if (program_number != blast_type_tblastn ||
          hit_params->options->longest_intron <= 0)
      {
         link_hsps(program_number, hsp_list, query_info, subject, sbp, 
                   hit_params, gapped_calculation);
         /* The HSP's may be in a different order than they were before, 
            but hsp contains the first one. */
      } else {
         /* Calculate individual HSP e-values first - they'll be needed to
            compare with sum e-values. */
         Blast_HSPListGetEvalues(program_number, query_info, 
                                    hsp_list, gapped_calculation, sbp);
         
         new_link_hsps(program_number, hsp_list, query_info, subject, sbp, 
                       hit_params);
      }
	}

	return 0;
}