link_hsps.c

ncbi源码

                  {
                     maxscore=sum;
                     best[1]=H;
                  }
               }
            }
            if(path_changed==0){
               /* No path was changed, use these max sums. */
               use_current_max=1;
            }
            else{
               /* If max path hasn't chaged, we can use it */
               /* Walk down best, give up if we find a removed item in path */
               use_current_max=1;
               if(!ignore_small_gaps){
                  for (H=best[0]; H!=NULL; H=H->hsp_link.link[0]) 
                     if (H->linked_to==-1000) {use_current_max=0; break;}
               }
               if(use_current_max)
                  for (H=best[1]; H!=NULL; H=H->hsp_link.link[1]) 
                     if (H->linked_to==-1000) {use_current_max=0; break;}
               
            }
         }

         /* reset helper_info */
         /* Inside this while loop, the linked list order never changes 
          * So here we initialize an array of commonly used info, 
          * and in this loop we access these arrays instead of the actual list 
          */
         if(!use_current_max){
            for (H=hp_start,H_index=1; H!=NULL; H=H->next,H_index++) {
               Int4 s_frame = H->hsp->subject.frame;
               Int4 s_off_t = H->s_offset_trim;
               Int4 q_off_t = H->q_offset_trim;
               lh_helper[H_index].ptr = H;
               lh_helper[H_index].q_off_trim = q_off_t;
               lh_helper[H_index].s_off_trim = s_off_t;
               for(i=0;i<BLAST_NUMBER_OF_ORDERING_METHODS;i++)
                  lh_helper[H_index].sum[i] = H->hsp_link.sum[i];
               max[SIGN(s_frame)+1]=
                  MAX(max[SIGN(s_frame)+1],H->hsp_link.sum[1]);
               lh_helper[H_index].maxsum1 =max[SIGN(s_frame)+1];					   
               
               /* set next_larger to link back to closest entry with a sum1 
                  larger than this */
               {
                  Int4 cur_sum=lh_helper[H_index].sum[1];
                  Int4 prev = H_index-1;
                  Int4 prev_sum = lh_helper[prev].sum[1];
                  while((cur_sum>=prev_sum) && (prev>0)){
                     prev=lh_helper[prev].next_larger;
                     prev_sum = lh_helper[prev].sum[1];
                  }
                  lh_helper[H_index].next_larger = prev;
               }
               H->linked_to = 0;
            }
            
            lh_helper[1].maxsum1 = -10000;
            
            /****** loop iter for index = 0  **************************/
            if(!ignore_small_gaps)
            {
               index=0;
               maxscore = -cutoff[index];
               H_index = 2;
               for (H=hp_start->next; H!=NULL; H=H->next,H_index++) 
               {
                  Int4 H_hsp_num=0;
                  Int4 H_hsp_sum=0;
                  double H_hsp_xsum=0.0;
                  LinkHSPStruct* H_hsp_link=NULL;
                  if (H->hsp->score > cutoff[index]) {
                     Int4 H_query_etrim = H->q_end_trim;
                     Int4 H_sub_etrim = H->s_end_trim;
                     Int4 H_q_et_gap = H_query_etrim+gap_size;
                     Int4 H_s_et_gap = H_sub_etrim+gap_size;
                     
                     /* We only walk down hits with the same frame sign */
                     /* for (H2=H->prev; H2!=NULL; H2=H2->prev,H2_index--) */
                     for (H2_index=H_index-1; H2_index>1; H2_index=H2_index-1) 
                     {
                        Int4 b1,b2,b4,b5;
                        Int4 q_off_t,s_off_t,sum;
                        
                        /* s_frame = lh_helper[H2_index].s_frame; */
                        q_off_t = lh_helper[H2_index].q_off_trim;
                        s_off_t = lh_helper[H2_index].s_off_trim;
                        
                        /* combine tests to reduce mispredicts -cfj */
                        b1 = q_off_t <= H_query_etrim;
                        b2 = s_off_t <= H_sub_etrim;
                        sum = lh_helper[H2_index].sum[index];
                        
                        
                        b4 = ( q_off_t > H_q_et_gap ) ;
                        b5 = ( s_off_t > H_s_et_gap ) ;
                        
                        /* list is sorted by q_off, so q_off should only increase.
                         * q_off_t can only differ from q_off by max_q_diff
                         * So once q_off_t is large enough (ie it exceeds limit 
                         * by max_q_diff), we can stop.  -cfj 
                         */
                        if(q_off_t > (H_q_et_gap+max_q_diff))
                           break;
                        
                        if (b1|b2|b5|b4) continue;
                        
                        if (sum>H_hsp_sum) 
                        {
                           H2=lh_helper[H2_index].ptr; 
                           H_hsp_num=H2->hsp_link.num[index];
                           H_hsp_sum=H2->hsp_link.sum[index];
                           H_hsp_xsum=H2->hsp_link.xsum[index];
                           H_hsp_link=H2;
                        }
                     } /* end for H2... */
                  }
                  { 
                     Int4 score=H->hsp->score;
                     double new_xsum = H_hsp_xsum + (score*(kbp[H->hsp->context]->Lambda));
                     Int4 new_sum = H_hsp_sum + (score - cutoff[index]);
                     
                     H->hsp_link.sum[index] = new_sum;
                     H->hsp_link.num[index] = H_hsp_num+1;
                     H->hsp_link.link[index] = H_hsp_link;
                     lh_helper[H_index].sum[index] = new_sum;
                     if (new_sum >= maxscore) 
                     {
                        maxscore=new_sum;
                        best[index]=H;
                     }
                     H->hsp_link.xsum[index] = new_xsum;
                     if(H_hsp_link)
                        ((LinkHSPStruct*)H_hsp_link)->linked_to++;
                  }
               } /* end for H=... */
            }
            /****** loop iter for index = 1  **************************/
            index=1;
            maxscore = -cutoff[index];
            H_index = 2;
            for (H=hp_start->next; H!=NULL; H=H->next,H_index++) 
            {
               Int4 H_hsp_num=0;
               Int4 H_hsp_sum=0;
               double H_hsp_xsum=0.0;
               LinkHSPStruct* H_hsp_link=NULL;
               
               H->hsp_link.changed=1;
               H2 = H->hsp_link.link[index];
               if ((!first_pass) && ((H2==0) || (H2->hsp_link.changed==0)))
               {
                  /* If the best choice last time has not been changed, then 
                     it is still the best choice, so no need to walk down list.
                  */
                  if(H2){
                     H_hsp_num=H2->hsp_link.num[index];
                     H_hsp_sum=H2->hsp_link.sum[index];
                     H_hsp_xsum=H2->hsp_link.xsum[index];
                  }
                  H_hsp_link=H2;
                  H->hsp_link.changed=0;
               } else if (H->hsp->score > cutoff[index]) {
                  Int4 H_query_etrim = H->q_end_trim;
                  Int4 H_sub_etrim = H->s_end_trim;
                  

                  /* Here we look at what was the best choice last time, if it's
                   * still around, and set this to the initial choice. By
                   * setting the best score to a (potentially) large value
                   * initially, we can reduce the number of hsps checked. 
                   */
                  
                  /* Currently we set the best score to a value just less than
                   * the real value. This is not really necessary, but doing
                   * this ensures that in the case of a tie, we make the same
                   * selection the original code did.
                   */
                  
                  if(!first_pass&&H2&&H2->linked_to>=0){
                     if(1){
                        /* We set this to less than the real value to keep the
                           original ordering in case of ties. */
                        H_hsp_sum=H2->hsp_link.sum[index]-1;
                     }else{
                        H_hsp_num=H2->hsp_link.num[index];
                        H_hsp_sum=H2->hsp_link.sum[index];
                        H_hsp_xsum=H2->hsp_link.xsum[index];
                        H_hsp_link=H2;
                     }
                  }
                  
                  /* We now only walk down hits with the same frame sign */
                  /* for (H2=H->prev; H2!=NULL; H2=H2->prev,H2_index--) */
                  for (H2_index=H_index-1; H2_index>1;)
                  {
                     Int4 b0,b1,b2;
                     Int4 q_off_t,s_off_t,sum,next_larger;
                     LinkHelpStruct * H2_helper=&lh_helper[H2_index];
                     sum = H2_helper->sum[index];
                     next_larger = H2_helper->next_larger;
                     
                     s_off_t = H2_helper->s_off_trim;
                     q_off_t = H2_helper->q_off_trim;
                     
                     b0 = sum <= H_hsp_sum;
                     
                     /* Compute the next H2_index */
                     H2_index--;
                     if(b0){	 /* If this sum is too small to beat H_hsp_sum, advance to a larger sum */
                        H2_index=next_larger;
                     }

                     /* combine tests to reduce mispredicts -cfj */
                     b1 = q_off_t <= H_query_etrim;
                     b2 = s_off_t <= H_sub_etrim;
                     
                     if(0) if(H2_helper->maxsum1<=H_hsp_sum)break;
                     
                     if (!(b0|b1|b2) )
                     {
                        H2 = H2_helper->ptr;
                        
                        H_hsp_num=H2->hsp_link.num[index];
                        H_hsp_sum=H2->hsp_link.sum[index];
                        H_hsp_xsum=H2->hsp_link.xsum[index];
                        H_hsp_link=H2;
                     }
                  } /* end for H2_index... */
               } /* end if(H->score>cuttof[]) */
               { 
                  Int4 score=H->hsp->score;
                  double new_xsum = 
                     H_hsp_xsum + (score*(kbp[H->hsp->context]->Lambda));
                  Int4 new_sum = H_hsp_sum + (score - cutoff[index]);
                  
                  H->hsp_link.sum[index] = new_sum;
                  H->hsp_link.num[index] = H_hsp_num+1;
                  H->hsp_link.link[index] = H_hsp_link;
                  lh_helper[H_index].sum[index] = new_sum;
                  lh_helper[H_index].maxsum1 = MAX(lh_helper[H_index-1].maxsum1, new_sum);
                  /* Update this entry's 'next_larger' field */
                  {
                     Int4 cur_sum=lh_helper[H_index].sum[1];
                     Int4 prev = H_index-1;
                     Int4 prev_sum = lh_helper[prev].sum[1];
                     while((cur_sum>=prev_sum) && (prev>0)){
                        prev=lh_helper[prev].next_larger;
                        prev_sum = lh_helper[prev].sum[1];
                     }
                     lh_helper[H_index].next_larger = prev;
                  }
                  
                  if (new_sum >= maxscore) 
                  {
                     maxscore=new_sum;
                     best[index]=H;
                  }
                  H->hsp_link.xsum[index] = new_xsum;
                  if(H_hsp_link)
                     ((LinkHSPStruct*)H_hsp_link)->linked_to++;
               }
            }
            path_changed=0;
            first_pass=0;
         }
         /********************************/
         if (!ignore_small_gaps)
         {
            /* Select the best ordering method.
               First we add back in the value cutoff[index] * the number
               of links, as this was subtracted out for purposes of the
               comparison above. */
            best[0]->hsp_link.sum[0] +=
               (best[0]->hsp_link.num[0])*cutoff[0];

            prob[0] = BLAST_SmallGapSumE(kbp[query_context],
                         gap_size,
                         best[0]->hsp_link.num[0], best[0]->hsp_link.xsum[0],
                         query_length, subject_length,
                         BLAST_GapDecayDivisor(gap_decay_rate,
                                              best[0]->hsp_link.num[0]) );

            /* Adjust the e-value because we are performing multiple tests */
            if( best[0]->hsp_link.num[0] > 1 ) {
              if( gap_prob == 0 || (prob[0] /= gap_prob) > INT4_MAX ) {
                prob[0] = INT4_MAX;
              }
            }

            prob[1] = BLAST_LargeGapSumE(kbp[query_context],
                         best[1]->hsp_link.num[1],
                         best[1]->hsp_link.xsum[1],
                         query_length, subject_length,
                         BLAST_GapDecayDivisor(gap_decay_rate,
                                              best[1]->hsp_link.num[1]));

            if( best[1]->hsp_link.num[1] > 1 ) {
              if( 1 - gap_prob == 0 || (prob[1] /= 1 - gap_prob) > INT4_MAX ) {
                prob[1] = INT4_MAX;
              }
            }
            ordering_method =
               prob[0]<=prob[1] ? BLAST_SMALL_GAPS : BLAST_LARGE_GAPS;
         }
         else
         {
            /* We only consider the case of big gaps. */
            best[1]->hsp_link.sum[1] +=
               (best[1]->hsp_link.num[1])*cutoff[1];

            prob[1] = BLAST_LargeGapSumE(kbp[query_context],
                         best[1]->hsp_link.num[1],
                         best[1]->hsp_link.xsum[1],
                         query_length, subject_length,
                         BLAST_GapDecayDivisor(gap_decay_rate,
                                              best[1]->hsp_link.num[1]));
            ordering_method = BLAST_LARGE_GAPS;
         }

         best[ordering_method]->start_of_chain = TRUE;
         
         /* AM: Support for query concatenation. */
         prob[ordering_method] *= 
            ((double)query_info->eff_searchsp_array[query_context] /
             ((double)subject_length*query_length));
         
         best[ordering_method]->hsp->evalue = prob[ordering_method];
         
         /* remove the links that have been ordered already. */
         if (best[ordering_method]->hsp_link.link[ordering_method])
            linked_set = TRUE;
         else
            linked_set = FALSE;

         if (best[ordering_method]->linked_to>0) path_changed=1;
         for (H=best[ordering_method]; H!=NULL;
              H=H->hsp_link.link[ordering_method]) 
         {
            if (H->linked_to>1) path_changed=1;
            H->linked_to=-1000;
            H->hsp_link.changed=1;
            /* record whether this is part of a linked set. */
            H->linked_set = linked_set;
            H->ordering_method = ordering_method;
            H->hsp->evalue = prob[ordering_method];
            if (H->next)
               (H->next)->prev=H->prev;
            if (H->prev)
               (H->prev)->next=H->next;
            number_of_hsps--;
         }
         
      } /* end while num_hsps... */