开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 生物技术 > 查看源码
link_hsps.c
资源名称：ncbi_cxx_win--Jul_08_2004.zip　[点击查看]
上传用户：yhdzpy8989
上传日期：2007-06-13
资源大小：13604k
文件大小：48k
源码类别：

生物技术

开发平台：

C/C++

link_hsps.c：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: link_hsps.c,v $
  4.  * PRODUCTION Revision 1000.3  2004/06/01 18:08:05  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.35
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /* $Id: link_hsps.c,v 1000.3 2004/06/01 18:08:05 gouriano Exp $
  11.  * ===========================================================================
  12.  *
  13.  *                            PUBLIC DOMAIN NOTICE
  14.  *               National Center for Biotechnology Information
  15.  *
  16.  *  This software/database is a "United States Government Work" under the
  17.  *  terms of the United States Copyright Act.  It was written as part of
  18.  *  the author's offical duties as a United States Government employee and
  19.  *  thus cannot be copyrighted.  This software/database is freely available
  20.  *  to the public for use. The National Library of Medicine and the U.S.
  21.  *  Government have not placed any restriction on its use or reproduction.
  22.  *
  23.  *  Although all reasonable efforts have been taken to ensure the accuracy
  24.  *  and reliability of the software and data, the NLM and the U.S.
  25.  *  Government do not and cannot warrant the performance or results that
  26.  *  may be obtained by using this software or data. The NLM and the U.S.
  27.  *  Government disclaim all warranties, express or implied, including
  28.  *  warranties of performance, merchantability or fitness for any particular
  29.  *  purpose.
  30.  *
  31.  *  Please cite the author in any work or product based on this material.
  32.  *
  33.  * ===========================================================================
  34.  *
  35.  * Author: Ilya Dondoshansky
  36.  *
  37.  */
  39. /** @file link_hsps.c
  40.  * Functions to link with use of sum statistics
  41.  */
  43. static char const rcsid[] = 
  44.     "$Id: link_hsps.c,v 1000.3 2004/06/01 18:08:05 gouriano Exp $";
  46. #include <algo/blast/core/link_hsps.h>
  47. #include <algo/blast/core/blast_util.h>
  50. /** Methods used to "order" the HSP's. */
  51. #define BLAST_NUMBER_OF_ORDERING_METHODS 2
  53. typedef enum LinkOrderingMethod {
  54.    BLAST_SMALL_GAPS = 0,
  55.    BLAST_LARGE_GAPS
  56. } LinkOrderingMethod;
  58. /* Forward declaration */
  59. struct LinkHSPStruct;
  61. /** The following structure is used in "link_hsps" to decide between
  62.  * two different "gapping" models.  Here link is used to hook up
  63.  * a chain of HSP's, num is the number of links, and sum is the sum score.
  64.  * Once the best gapping model has been found, this information is
  65.  * transferred up to the LinkHSPStruct.  This structure should not be
  66.  * used outside of the function link_hsps.
  67. */
  68. typedef struct BlastHSPLink {
  69.    struct LinkHSPStruct* link[BLAST_NUMBER_OF_ORDERING_METHODS]; /**< Best 
  70.                                                choice of HSP to link with */
  71.    Int2 num[BLAST_NUMBER_OF_ORDERING_METHODS]; /**< number of HSP in the
  72.                                                   ordering. */
  73.    Int4 sum[BLAST_NUMBER_OF_ORDERING_METHODS]; /**< Sum-Score of HSP. */
  74.    double xsum[BLAST_NUMBER_OF_ORDERING_METHODS]; /**< Sum-Score of HSP,
  75.                                      multiplied by the appropriate Lambda. */
  76.    Int4 changed; /**< Has the link been changed since previous access? */
  77. } BlastHSPLink;
  79. /** Structure containing all information internal to the process of linking
  80.  * HSPs.
  81.  */
  82. typedef struct LinkHSPStruct {
  83.    BlastHSP* hsp;      /**< Specific HSP this structure corresponds to */
  84.    struct LinkHSPStruct* prev;         /**< Previous HSP in a set, if any */
  85.    struct LinkHSPStruct* next;         /**< Next HSP in a set, if any */ 
  86.    BlastHSPLink  hsp_link; /**< Auxiliary structure for keeping track of sum
  87.                               scores, etc. */
  88.    Boolean linked_set;     /**< Is this HSp part of a linked set? */
  89.    Boolean start_of_chain; /**< If TRUE, this HSP starts a chain along the
  90.                               "link" pointer. */
  91.    Int4 linked_to;         /**< Where this HSP is linked to? */
  92.    Int4 sumscore;          /**< Sumscore of a set of "linked" HSP's. */
  93.    Int2 ordering_method;   /**< Which method (max or no max for gaps) was 
  94.                               used for linking HSPs? */
  95.    Int4 q_offset_trim;     /**< Start of trimmed hsp in query */
  96.    Int4 q_end_trim;        /**< End of trimmed HSP in query */
  97.    Int4 s_offset_trim;     /**< Start of trimmed hsp in subject */
  98.    Int4 s_end_trim;        /**< End of trimmed HSP in subject */
  99. } LinkHSPStruct;
  103. #define WINDOW_SIZE 20
  104. static double 
  105. SumHSPEvalue(Uint1 program_number, BlastScoreBlk* sbp, 
  106.    BlastQueryInfo* query_info, BLAST_SequenceBlk* subject, 
  107.    const BlastHitSavingParameters* hit_params, 
  108.    LinkHSPStruct* head_hsp, LinkHSPStruct* hsp, Int4* sumscore)
  109. {
  110.    double gap_prob, gap_decay_rate, sum_evalue, score_prime;
  111.    Int2 num;
  112.    Int4 subject_eff_length, query_eff_length, length_adjustment;
  113.    Int4 context = head_hsp->hsp->context;
  114.    double eff_searchsp;
  116.    gap_prob = hit_params->gap_prob;
  117.    gap_decay_rate = hit_params->gap_decay_rate;
  119.    num = head_hsp->hsp->num + hsp->hsp->num;
  121.    length_adjustment = query_info->length_adjustments[context];
  123.    subject_eff_length = 
  124.       MAX((subject->length - length_adjustment), 1);
  125.    if (program_number == blast_type_tblastn) 
  126.       subject_eff_length /= 3;
  127.    subject_eff_length = MAX(subject_eff_length, 1);
  129.    query_eff_length = 
  130.       MAX(BLAST_GetQueryLength(query_info, context) - length_adjustment, 1);
  131.    
  132.    *sumscore = MAX(hsp->hsp->score, hsp->sumscore) + 
  133.       MAX(head_hsp->hsp->score, head_hsp->sumscore);
  134.    score_prime = *sumscore * sbp->kbp_gap[context]->Lambda;
  136.    sum_evalue =
  137.        BLAST_UnevenGapSumE(sbp->kbp_gap[context], 2*WINDOW_SIZE,
  138.                            hit_params->options->longest_intron + WINDOW_SIZE,
  139.                            num, score_prime,
  140.                            query_eff_length, subject_eff_length,
  141.                            BLAST_GapDecayDivisor(gap_decay_rate, num));
  143.    eff_searchsp = ((double) subject_eff_length) * query_eff_length;
  144.    
  145.    sum_evalue *= 
  146.       ((double) query_info->eff_searchsp_array[context]) / eff_searchsp;
  148.    return sum_evalue;
  149. }
  151. /** Sort the HSP's by starting position of the query.  Called by qsort.  
  152.  * The first function sorts in forward, the second in reverse order.
  153. */
  154. static int
  155. fwd_compare_hsps(const void* v1, const void* v2)
  156. {
  157. BlastHSP* h1,* h2;
  158. LinkHSPStruct** hp1,** hp2;
  160. hp1 = (LinkHSPStruct**) v1;
  161. hp2 = (LinkHSPStruct**) v2;
  162. h1 = (*hp1)->hsp;
  163. h2 = (*hp2)->hsp;
  165. if (h1->context < h2->context)
  166.       return -1;
  167.    else if (h1->context > h2->context)
  168.       return 1;
  170. if (h1->query.offset < h2->query.offset) 
  171. return -1;
  172. if (h1->query.offset > h2->query.offset) 
  173. return 1;
  174. /* Necessary in case both HSP's have the same query offset. */
  175. if (h1->subject.offset < h2->subject.offset) 
  176. return -1;
  177. if (h1->subject.offset > h2->subject.offset) 
  178. return 1;
  180. return 0;
  181. }
  183. /** Sort the HSP's by starting position of the query.  Called by qsort.  
  184.  * The first function sorts in forward, the second in reverse order.
  185. */
  186. static int
  187. fwd_compare_hsps_transl(const void* v1, const void* v2)
  188. {
  189. BlastHSP* h1,* h2;
  190. LinkHSPStruct** hp1,** hp2;
  191.    Int4 context1, context2;
  193. hp1 = (LinkHSPStruct**) v1;
  194. hp2 = (LinkHSPStruct**) v2;
  195. h1 = (*hp1)->hsp;
  196. h2 = (*hp2)->hsp;
  198.    context1 = h1->context/3;
  199.    context2 = h2->context/3;
  201.    if (context1 < context2)
  202.       return 1;
  203.    else if (context1 > context2)
  204.       return -1;
  206. if (h1->query.offset < h2->query.offset) 
  207. return -1;
  208. if (h1->query.offset > h2->query.offset) 
  209. return 1;
  210. /* Necessary in case both HSP's have the same query offset. */
  211. if (h1->subject.offset < h2->subject.offset) 
  212. return -1;
  213. if (h1->subject.offset > h2->subject.offset) 
  214. return 1;
  216. return 0;
  217. }
  219. /* Comparison function based on end position in the query */
  220. static int
  221. end_compare_hsps(const void* v1, const void* v2)
  222. {
  223. BlastHSP* h1,* h2;
  224. LinkHSPStruct** hp1,** hp2;
  226. hp1 = (LinkHSPStruct**) v1;
  227. hp2 = (LinkHSPStruct**) v2;
  228. h1 = (*hp1)->hsp;
  229. h2 = (*hp2)->hsp;
  231.    if (h1->context < h2->context)
  232.       return 1;
  233.    else if (h1->context > h2->context)
  234.       return -1;
  236. if (h1->query.end < h2->query.end) 
  237. return -1;
  238. if (h1->query.end > h2->query.end) 
  239. return 1;
  240. /* Necessary in case both HSP's have the same query end. */
  241. if (h1->subject.end < h2->subject.end) 
  242. return -1;
  243. if (h1->subject.end > h2->subject.end) 
  244. return 1;
  246. return 0;
  247. }
  249. static int
  250. sumscore_compare_hsps(const void* v1, const void* v2)
  251. {
  252. LinkHSPStruct* h1,* h2;
  253. LinkHSPStruct** hp1,** hp2;
  254.    Int4 score1, score2;
  256. hp1 = (LinkHSPStruct**) v1;
  257. hp2 = (LinkHSPStruct**) v2;
  258. h1 = *hp1;
  259. h2 = *hp2;
  261. if (h1 == NULL || h2 == NULL)
  262. return 0;
  264.    score1 = MAX(h1->sumscore, h1->hsp->score);
  265.    score2 = MAX(h2->sumscore, h2->hsp->score);
  267. if (score1 < score2) 
  268. return 1;
  269. if (score1 > score2)
  270. return -1;
  272. return 0;
  273. }
  276. /* The following function should be used only for new tblastn. 
  277.    Current implementation does not allow its use in two sequences BLAST */
  279. #define MAX_SPLICE_DIST 5
  280. static Boolean
  281. FindSpliceJunction(Uint1* subject_seq, BlastHSP* hsp1, BlastHSP* hsp2)
  282. {
  283.    Boolean found = FALSE;
  284.    Int4 overlap, length, i;
  285.    Uint1* nt_seq;
  286.    Uint1 g = 4, t = 8, a = 1; /* ncbi4na values for G, T, A respectively */
  288.    if (!subject_seq)
  289.       return FALSE;
  291.    overlap = hsp1->query.end - hsp2->query.offset;
  293.    if (overlap >= 0) {
  294.       length = 3*overlap + 2;
  295.       nt_seq = &subject_seq[hsp1->subject.end - 3*overlap];
  296.    } else {
  297.       length = MAX_SPLICE_DIST;
  298.       nt_seq = &subject_seq[hsp1->subject.end];
  299.    }
  300.    
  301.    for (i=0; i<length-1; i++) {
  302.       if (nt_seq[i] == g && nt_seq[i+1] == t) {
  303.          found = TRUE;
  304.          break;
  305.       }
  306.    }
  307.       
  308.    if (!found) 
  309.       return FALSE;
  310.    else
  311.       found = FALSE;
  313.    if (overlap >= 0) 
  314.       nt_seq = &subject_seq[hsp2->subject.offset - 2];
  315.    else 
  316.       nt_seq = &subject_seq[hsp2->subject.offset - MAX_SPLICE_DIST];
  318.    for (i=0; i<length-1; i++) {
  319.       if (nt_seq[i] == a && nt_seq[i+1] == g) {
  320.          found = TRUE;
  321.          break;
  322.       }
  323.    }
  324.    return found;
  325. }
  327. /* Find an HSP with offset closest, but not smaller/larger than a given one.
  328.  */
  329. static Int4 hsp_binary_search(LinkHSPStruct** hsp_array, Int4 size, 
  330.                               Int4 offset, Boolean right)
  331. {
  332.    Int4 index, begin, end, coord;
  333.    
  334.    begin = 0;
  335.    end = size;
  336.    while (begin < end) {
  337.       index = (begin + end) / 2;
  338.       if (right) 
  339.          coord = hsp_array[index]->hsp->query.offset;
  340.       else
  341.          coord = hsp_array[index]->hsp->query.end;
  342.       if (coord >= offset) 
  343.          end = index;
  344.       else
  345.          begin = index + 1;
  346.    }
  348.    return end;
  349. }
  351. static int
  352. rev_compare_hsps(const void *v1, const void *v2)
  354. {
  355. BlastHSP* h1,* h2;
  356. LinkHSPStruct** hp1,** hp2;
  358. hp1 = (LinkHSPStruct**) v1;
  359. hp2 = (LinkHSPStruct**) v2;
  360. h1 = (*hp1)->hsp;
  361. h2 = (*hp2)->hsp;
  363.    if (h1->context > h2->context)
  364.       return 1;
  365.    else if (h1->context < h2->context)
  366.       return -1;
  368. if (h1->query.offset < h2->query.offset) 
  369. return  1;
  370. if (h1->query.offset > h2->query.offset) 
  371. return -1;
  372. return 0;
  373. }
  375. static int
  376. rev_compare_hsps_transl(const void *v1, const void *v2)
  378. {
  379. BlastHSP* h1,* h2;
  380. LinkHSPStruct** hp1,** hp2;
  381.    Int4 context1, context2;
  383. hp1 = (LinkHSPStruct**) v1;
  384. hp2 = (LinkHSPStruct**) v2;
  385. h1 = (*hp1)->hsp;
  386. h2 = (*hp2)->hsp;
  388.    context1 = h1->context/3;
  389.    context2 = h2->context/3;
  391.    if (context1 > context2)
  392.       return 1;
  393.    else if (context1 < context2)
  394.       return -1;
  396. if (h1->query.offset < h2->query.offset) 
  397. return  1;
  398. if (h1->query.offset > h2->query.offset) 
  399. return -1;
  400. return 0;
  401. }
  403. static int
  404. rev_compare_hsps_tbn(const void *v1, const void *v2)
  406. {
  407. BlastHSP* h1,* h2;
  408. LinkHSPStruct** hp1,** hp2;
  410. hp1 = (LinkHSPStruct**) v1;
  411. hp2 = (LinkHSPStruct**) v2;
  412. h1 = (*hp1)->hsp;
  413. h2 = (*hp2)->hsp;
  415.    if (h1->context > h2->context)
  416.       return 1;
  417.    else if (h1->context < h2->context)
  418.       return -1;
  420. if (SIGN(h1->subject.frame) != SIGN(h2->subject.frame))
  421. {
  422. if (h1->subject.frame > h2->subject.frame)
  423. return 1;
  424. else
  425. return -1;
  426. }
  428. if (h1->query.offset < h2->query.offset) 
  429. return  1;
  430. if (h1->query.offset > h2->query.offset) 
  431. return -1;
  432. return 0;
  433. }
  435. static int
  436. rev_compare_hsps_tbx(const void *v1, const void *v2)
  438. {
  439. BlastHSP* h1,* h2;
  440. LinkHSPStruct** hp1,** hp2;
  441.    Int4 context1, context2;
  443. hp1 = (LinkHSPStruct**) v1;
  444. hp2 = (LinkHSPStruct**) v2;
  445. h1 = (*hp1)->hsp;
  446. h2 = (*hp2)->hsp;
  448.    context1 = h1->context/3;
  449.    context2 = h2->context/3;
  451.    if (context1 > context2)
  452.       return 1;
  453.    else if (context1 < context2)
  454.       return -1;
  455.    
  456. if (SIGN(h1->subject.frame) != SIGN(h2->subject.frame))
  457. {
  458. if (h1->subject.frame > h2->subject.frame)
  459. return 1;
  460. else
  461. return -1;
  462. }
  464. if (h1->query.offset < h2->query.offset) 
  465. return  1;
  466. if (h1->query.offset > h2->query.offset) 
  467. return -1;
  468. return 0;
  469. }
  471. /** The helper array contains the info used frequently in the inner 
  472.  * for loops of the HSP linking algorithm.
  473.  * One array of helpers will be allocated for each thread.
  474.  */
  475. typedef struct LinkHelpStruct {
  476.   LinkHSPStruct* ptr;
  477.   Int4 q_off_trim;
  478.   Int4 s_off_trim;
  479.   Int4 sum[BLAST_NUMBER_OF_ORDERING_METHODS];
  480.   Int4 maxsum1;
  481.   Int4 next_larger;
  482. } LinkHelpStruct;
  484. static LinkHSPStruct* LinkHSPStructReset(LinkHSPStruct* lhsp)
  485. {
  486.    BlastHSP* hsp;
  488.    if (!lhsp) {
  489.       lhsp = (LinkHSPStruct*) calloc(1, sizeof(LinkHSPStruct));
  490.       lhsp->hsp = (BlastHSP*) calloc(1, sizeof(BlastHSP));
  491.    } else {
  492.       if (!lhsp->hsp) {
  493.          hsp = (BlastHSP*) calloc(1, sizeof(BlastHSP));
  494.       } else {
  495.          hsp = lhsp->hsp;
  496.          memset(hsp, 0, sizeof(BlastHSP));
  497.       }
  498.       memset(lhsp, 0, sizeof(LinkHSPStruct));
  499.       lhsp->hsp = hsp;
  500.    }
  501.    return lhsp;
  502. }
  504. static Int2
  505. link_hsps(Uint1 program_number, BlastHSPList* hsp_list, 
  506.    BlastQueryInfo* query_info, BLAST_SequenceBlk* subject,
  507.    BlastScoreBlk* sbp, const BlastHitSavingParameters* hit_params,
  508.    Boolean gapped_calculation)
  509. {
  510. LinkHSPStruct* H,* H2,* best[2],* first_hsp,* last_hsp,** hp_frame_start;
  511. LinkHSPStruct* hp_start = NULL;
  512.    BlastHSP* hsp;
  513.    BlastHSP** hsp_array;
  514. Blast_KarlinBlk** kbp;
  515. Int4 maxscore, cutoff[2];
  516. Boolean linked_set, ignore_small_gaps;
  517. double gap_decay_rate, gap_prob, prob[2];
  518. Int4 index, index1, num_links, frame_index;
  519.    LinkOrderingMethod ordering_method;
  520.    Int4 num_query_frames, num_subject_frames;
  521. Int4 *hp_frame_number;
  522. Int4 gap_size, number_of_hsps, total_number_of_hsps;
  523.    Int4 query_length, subject_length, length_adjustment;
  524. LinkHSPStruct* link;
  525. Int4 H2_index,H_index;
  526. Int4 i;
  527. Int4 max_q_diff;
  528.   Int4 path_changed;  /* will be set if an element is removed that may change an existing path */
  529.   Int4 first_pass, use_current_max; 
  530. LinkHelpStruct *lh_helper=0;
  531.    Int4 lh_helper_size;
  532. Int4 query_context; /* AM: to support query concatenation. */
  533.    Boolean translated_query;
  534.    LinkHSPStruct** link_hsp_array;
  536. if (hsp_list == NULL)
  537. return -1;
  539.    hsp_array = hsp_list->hsp_array;
  541.    lh_helper_size = MAX(1024,hsp_list->hspcnt+5);
  542.    lh_helper = (LinkHelpStruct *) 
  543.       calloc(lh_helper_size, sizeof(LinkHelpStruct));
  545. if (gapped_calculation) 
  546. kbp = sbp->kbp_gap;
  547. else
  548. kbp = sbp->kbp;
  550. total_number_of_hsps = hsp_list->hspcnt;
  552. number_of_hsps = total_number_of_hsps;
  553. gap_size = hit_params->gap_size;
  554. gap_prob = hit_params->gap_prob;
  555. gap_decay_rate = hit_params->gap_decay_rate;
  557. translated_query = (program_number == blast_type_blastx || 
  558.                        program_number == blast_type_tblastx);
  559.    if (program_number == blast_type_tblastn ||
  560.        program_number == blast_type_tblastx)
  561.       num_subject_frames = NUM_STRANDS;
  562.    else
  563.       num_subject_frames = 1;
  565.    link_hsp_array = 
  566.       (LinkHSPStruct**) malloc(total_number_of_hsps*sizeof(LinkHSPStruct*));
  567.    for (index = 0; index < total_number_of_hsps; ++index) {
  568.       link_hsp_array[index] = (LinkHSPStruct*) calloc(1, sizeof(LinkHSPStruct));
  569.       link_hsp_array[index]->hsp = hsp_array[index];
  570.    }
  572.    /* Sort by (reverse) position. */
  573.    if (translated_query) {
  574.       qsort(link_hsp_array,total_number_of_hsps,sizeof(LinkHSPStruct*), 
  575.             rev_compare_hsps_tbx);
  576.    } else {
  577.       qsort(link_hsp_array,total_number_of_hsps,sizeof(LinkHSPStruct*), 
  578.             rev_compare_hsps_tbn);
  579.    }
  581. cutoff[0] = hit_params->cutoff_small_gap;
  582. cutoff[1] = hit_params->cutoff_big_gap;
  583. ignore_small_gaps = hit_params->ignore_small_gaps;
  585. if (translated_query)
  586. num_query_frames = NUM_STRANDS*query_info->num_queries;
  587. else
  588. num_query_frames = query_info->num_queries;
  589.    
  590.    hp_frame_start = 
  591.        calloc(num_subject_frames*num_query_frames, sizeof(LinkHSPStruct*));
  592.    hp_frame_number = calloc(num_subject_frames*num_query_frames, sizeof(Int4));
  594. /* hook up the HSP's */
  595. hp_frame_start[0] = link_hsp_array[0];
  597. /* Put entries with different frame parity into separate 'query_frame's. -cfj */
  598. {
  599.   Int4 cur_frame=0;
  600.      for (index=0;index<number_of_hsps;index++) 
  601.      {
  602.         H=link_hsp_array[index];
  603.         H->start_of_chain = FALSE;
  604.         hp_frame_number[cur_frame]++;
  606.         H->prev= index ? link_hsp_array[index-1] : NULL;
  607.         H->next= index<(number_of_hsps-1) ? link_hsp_array[index+1] : NULL;
  608.         if (H->prev != NULL && 
  609.             ((H->hsp->context/3) != (H->prev->hsp->context/3) ||
  610.              (SIGN(H->hsp->subject.frame) != SIGN(H->prev->hsp->subject.frame))))
  611.         { /* If frame switches, then start new list. */
  612.            hp_frame_number[cur_frame]--;
  613.            hp_frame_number[++cur_frame]++;
  614.            hp_frame_start[cur_frame] = H;
  615.            H->prev->next = NULL;
  616.            H->prev = NULL;
  617.         }
  618.      }
  619.      num_query_frames = cur_frame+1;
  620. }
  622. /* max_q_diff is the maximum amount q.offset can differ from q.offset_trim */
  623. /* This is used to break out of H2 loop early */
  624.    for (index=0;index<number_of_hsps;index++) 
  625.    {
  626.       H = link_hsp_array[index];
  627. hsp = H->hsp;
  628. H->q_offset_trim = hsp->query.offset + MIN(((hsp->query.length)/4), 5);
  629. H->q_end_trim = hsp->query.end - MIN(((hsp->query.length)/4), 5);
  630. H->s_offset_trim = 
  631.          hsp->subject.offset + MIN(((hsp->subject.length)/4), 5);
  632. H->s_end_trim = hsp->subject.end - MIN(((hsp->subject.length)/4), 5);
  633.    }     
  634.    max_q_diff = 5;
  636. for (frame_index=0; frame_index<num_query_frames; frame_index++)
  637. {
  638.       hp_start = LinkHSPStructReset(hp_start);
  639.       hp_start->next = hp_frame_start[frame_index];
  640.       hp_frame_start[frame_index]->prev = hp_start;
  641.       number_of_hsps = hp_frame_number[frame_index];
  642.       query_context = hp_start->next->hsp->context;
  643.       length_adjustment = query_info->length_adjustments[query_context];
  644.       query_length = BLAST_GetQueryLength(query_info, query_context);
  645.       query_length = MAX(query_length - length_adjustment, 1);
  646.       subject_length = subject->length; /* in nucleotides even for tblast[nx] */
  647.       /* If subject is translated, length adjustment is given in nucleotide
  648.          scale. */
  649.       if (program_number == blast_type_tblastn || 
  650.           program_number == blast_type_tblastx)
  651.       {
  652.          length_adjustment /= CODON_LENGTH;
  653.          subject_length /= CODON_LENGTH;
  654.       }
  655.       subject_length = MAX(subject_length - length_adjustment, 1);
  657.       lh_helper[0].ptr = hp_start;
  658.       lh_helper[0].q_off_trim = 0;
  659.       lh_helper[0].s_off_trim = 0;
  660.       lh_helper[0].maxsum1  = -10000;
  661.       lh_helper[0].next_larger  = 0;
  662.       
  663.       /* lh_helper[0]  = empty     = additional end marker
  664.        * lh_helper[1]  = hsp_start = empty entry used in original code
  665.        * lh_helper[2]  = hsp_array->next = hsp_array[0]
  666.        * lh_helper[i]  = ... = hsp_array[i-2] (for i>=2) 
  667.        */
  668.       first_pass=1;    /* do full search */
  669.       path_changed=1;
  670.       for (H=hp_start->next; H!=NULL; H=H->next) 
  671.          H->hsp_link.changed=1;
  673.       while (number_of_hsps > 0)
  674.       {
  675.          Int4 max[3];
  676.          max[0]=max[1]=max[2]=-10000;
  677.          /* Initialize the 'best' parameter */
  678.          best[0] = best[1] = NULL;
  679.          
  680.          
  681.          /* See if we can avoid recomputing all scores:
  682.           *  - Find the max paths (based on old scores). 
  683.           *  - If no paths were changed by removal of nodes (ie research==0) 
  684.           *    then these max paths are still the best.
  685.           *  - else if these max paths were unchanged, then they are still the best.
  686.           */
  687.          use_current_max=0;
  688.          if (!first_pass){
  689.             Int4 max0,max1;
  690.             /* Find the current max sums */
  691.             if(!ignore_small_gaps){
  692.                max0 = -cutoff[0];
  693.                max1 = -cutoff[1];
  694.                for (H=hp_start->next; H!=NULL; H=H->next) {
  695.                   Int4 sum0=H->hsp_link.sum[0];
  696.                   Int4 sum1=H->hsp_link.sum[1];
  697.                   if(sum0>=max0)
  698.                   {
  699.                      max0=sum0;
  700.                      best[0]=H;
  701.                   }
  702.                   if(sum1>=max1)
  703.                   {
  704.                      max1=sum1;
  705.                      best[1]=H;
  706.                   }
  707.                }
  708.             } else {
  709.                maxscore = -cutoff[1];
  710.                for (H=hp_start->next; H!=NULL; H=H->next) {
  711.                   Int4  sum=H->hsp_link.sum[1];
  712.                   if(sum>=maxscore)
  713.                   {
  714.                      maxscore=sum;
  715.                      best[1]=H;
  716.                   }
  717.                }
  718.             }
  719.             if(path_changed==0){
  720.                /* No path was changed, use these max sums. */
  721.                use_current_max=1;
  722.             }
  723.             else{
  724.                /* If max path hasn't chaged, we can use it */
  725.                /* Walk down best, give up if we find a removed item in path */
  726.                use_current_max=1;
  727.                if(!ignore_small_gaps){
  728.                   for (H=best[0]; H!=NULL; H=H->hsp_link.link[0]) 
  729.                      if (H->linked_to==-1000) {use_current_max=0; break;}
  730.                }
  731.                if(use_current_max)
  732.                   for (H=best[1]; H!=NULL; H=H->hsp_link.link[1]) 
  733.                      if (H->linked_to==-1000) {use_current_max=0; break;}
  734.                
  735.             }
  736.          }
  738.          /* reset helper_info */
  739.          /* Inside this while loop, the linked list order never changes 
  740.           * So here we initialize an array of commonly used info, 
  741.           * and in this loop we access these arrays instead of the actual list 
  742.           */
  743.          if(!use_current_max){
  744.             for (H=hp_start,H_index=1; H!=NULL; H=H->next,H_index++) {
  745.                Int4 s_frame = H->hsp->subject.frame;
  746.                Int4 s_off_t = H->s_offset_trim;
  747.                Int4 q_off_t = H->q_offset_trim;
  748.                lh_helper[H_index].ptr = H;
  749.                lh_helper[H_index].q_off_trim = q_off_t;
  750.                lh_helper[H_index].s_off_trim = s_off_t;
  751.                for(i=0;i<BLAST_NUMBER_OF_ORDERING_METHODS;i++)
  752.                   lh_helper[H_index].sum[i] = H->hsp_link.sum[i];
  753.                max[SIGN(s_frame)+1]=
  754.                   MAX(max[SIGN(s_frame)+1],H->hsp_link.sum[1]);
  755.                lh_helper[H_index].maxsum1 =max[SIGN(s_frame)+1];    
  756.                
  757.                /* set next_larger to link back to closest entry with a sum1 
  758.                   larger than this */
  759.                {
  760.                   Int4 cur_sum=lh_helper[H_index].sum[1];
  761.                   Int4 prev = H_index-1;
  762.                   Int4 prev_sum = lh_helper[prev].sum[1];
  763.                   while((cur_sum>=prev_sum) && (prev>0)){
  764.                      prev=lh_helper[prev].next_larger;
  765.                      prev_sum = lh_helper[prev].sum[1];
  766.                   }
  767.                   lh_helper[H_index].next_larger = prev;
  768.                }
  769.                H->linked_to = 0;
  770.             }
  771.             
  772.             lh_helper[1].maxsum1 = -10000;
  773.             
  774.             /****** loop iter for index = 0  **************************/
  775.             if(!ignore_small_gaps)
  776.             {
  777.                index=0;
  778.                maxscore = -cutoff[index];
  779.                H_index = 2;
  780.                for (H=hp_start->next; H!=NULL; H=H->next,H_index++) 
  781.                {
  782.                   Int4 H_hsp_num=0;
  783.                   Int4 H_hsp_sum=0;
  784.                   double H_hsp_xsum=0.0;
  785.                   LinkHSPStruct* H_hsp_link=NULL;
  786.                   if (H->hsp->score > cutoff[index]) {
  787.                      Int4 H_query_etrim = H->q_end_trim;
  788.                      Int4 H_sub_etrim = H->s_end_trim;
  789.                      Int4 H_q_et_gap = H_query_etrim+gap_size;
  790.                      Int4 H_s_et_gap = H_sub_etrim+gap_size;
  791.                      
  792.                      /* We only walk down hits with the same frame sign */
  793.                      /* for (H2=H->prev; H2!=NULL; H2=H2->prev,H2_index--) */
  794.                      for (H2_index=H_index-1; H2_index>1; H2_index=H2_index-1) 
  795.                      {
  796.                         Int4 b1,b2,b4,b5;
  797.                         Int4 q_off_t,s_off_t,sum;
  798.                         
  799.                         /* s_frame = lh_helper[H2_index].s_frame; */
  800.                         q_off_t = lh_helper[H2_index].q_off_trim;
  801.                         s_off_t = lh_helper[H2_index].s_off_trim;
  802.                         
  803.                         /* combine tests to reduce mispredicts -cfj */
  804.                         b1 = q_off_t <= H_query_etrim;
  805.                         b2 = s_off_t <= H_sub_etrim;
  806.                         sum = lh_helper[H2_index].sum[index];
  807.                         
  808.                         
  809.                         b4 = ( q_off_t > H_q_et_gap ) ;
  810.                         b5 = ( s_off_t > H_s_et_gap ) ;
  811.                         
  812.                         /* list is sorted by q_off, so q_off should only increase.
  813.                          * q_off_t can only differ from q_off by max_q_diff
  814.                          * So once q_off_t is large enough (ie it exceeds limit 
  815.                          * by max_q_diff), we can stop.  -cfj 
  816.                          */
  817.                         if(q_off_t > (H_q_et_gap+max_q_diff))
  818.                            break;
  819.                         
  820.                         if (b1|b2|b5|b4) continue;
  821.                         
  822.                         if (sum>H_hsp_sum) 
  823.                         {
  824.                            H2=lh_helper[H2_index].ptr; 
  825.                            H_hsp_num=H2->hsp_link.num[index];
  826.                            H_hsp_sum=H2->hsp_link.sum[index];
  827.                            H_hsp_xsum=H2->hsp_link.xsum[index];
  828.                            H_hsp_link=H2;
  829.                         }
  830.                      } /* end for H2... */
  831.                   }
  832.                   { 
  833.                      Int4 score=H->hsp->score;
  834.                      double new_xsum = H_hsp_xsum + (score*(kbp[H->hsp->context]->Lambda));
  835.                      Int4 new_sum = H_hsp_sum + (score - cutoff[index]);
  836.                      
  837.                      H->hsp_link.sum[index] = new_sum;
  838.                      H->hsp_link.num[index] = H_hsp_num+1;
  839.                      H->hsp_link.link[index] = H_hsp_link;
  840.                      lh_helper[H_index].sum[index] = new_sum;
  841.                      if (new_sum >= maxscore) 
  842.                      {
  843.                         maxscore=new_sum;
  844.                         best[index]=H;
  845.                      }
  846.                      H->hsp_link.xsum[index] = new_xsum;
  847.                      if(H_hsp_link)
  848.                         ((LinkHSPStruct*)H_hsp_link)->linked_to++;
  849.                   }
  850.                } /* end for H=... */
  851.             }
  852.             /****** loop iter for index = 1  **************************/
  853.             index=1;
  854.             maxscore = -cutoff[index];
  855.             H_index = 2;
  856.             for (H=hp_start->next; H!=NULL; H=H->next,H_index++) 
  857.             {
  858.                Int4 H_hsp_num=0;
  859.                Int4 H_hsp_sum=0;
  860.                double H_hsp_xsum=0.0;
  861.                LinkHSPStruct* H_hsp_link=NULL;
  862.                
  863.                H->hsp_link.changed=1;
  864.                H2 = H->hsp_link.link[index];
  865.                if ((!first_pass) && ((H2==0) || (H2->hsp_link.changed==0)))
  866.                {
  867.                   /* If the best choice last time has not been changed, then 
  868.                      it is still the best choice, so no need to walk down list.
  869.                   */
  870.                   if(H2){
  871.                      H_hsp_num=H2->hsp_link.num[index];
  872.                      H_hsp_sum=H2->hsp_link.sum[index];
  873.                      H_hsp_xsum=H2->hsp_link.xsum[index];
  874.                   }
  875.                   H_hsp_link=H2;
  876.                   H->hsp_link.changed=0;
  877.                } else if (H->hsp->score > cutoff[index]) {
  878.                   Int4 H_query_etrim = H->q_end_trim;
  879.                   Int4 H_sub_etrim = H->s_end_trim;
  880.                   
  882.                   /* Here we look at what was the best choice last time, if it's
  883.                    * still around, and set this to the initial choice. By
  884.                    * setting the best score to a (potentially) large value
  885.                    * initially, we can reduce the number of hsps checked. 
  886.                    */
  887.                   
  888.                   /* Currently we set the best score to a value just less than
  889.                    * the real value. This is not really necessary, but doing
  890.                    * this ensures that in the case of a tie, we make the same
  891.                    * selection the original code did.
  892.                    */
  893.                   
  894.                   if(!first_pass&&H2&&H2->linked_to>=0){
  895.                      if(1){
  896.                         /* We set this to less than the real value to keep the
  897.                            original ordering in case of ties. */
  898.                         H_hsp_sum=H2->hsp_link.sum[index]-1;
  899.                      }else{
  900.                         H_hsp_num=H2->hsp_link.num[index];
  901.                         H_hsp_sum=H2->hsp_link.sum[index];
  902.                         H_hsp_xsum=H2->hsp_link.xsum[index];
  903.                         H_hsp_link=H2;
  904.                      }
  905.                   }
  906.                   
  907.                   /* We now only walk down hits with the same frame sign */
  908.                   /* for (H2=H->prev; H2!=NULL; H2=H2->prev,H2_index--) */
  909.                   for (H2_index=H_index-1; H2_index>1;)
  910.                   {
  911.                      Int4 b0,b1,b2;
  912.                      Int4 q_off_t,s_off_t,sum,next_larger;
  913.                      LinkHelpStruct * H2_helper=&lh_helper[H2_index];
  914.                      sum = H2_helper->sum[index];
  915.                      next_larger = H2_helper->next_larger;
  916.                      
  917.                      s_off_t = H2_helper->s_off_trim;
  918.                      q_off_t = H2_helper->q_off_trim;
  919.                      
  920.                      b0 = sum <= H_hsp_sum;
  921.                      
  922.                      /* Compute the next H2_index */
  923.                      H2_index--;
  924.                      if(b0){  /* If this sum is too small to beat H_hsp_sum, advance to a larger sum */
  925.                         H2_index=next_larger;
  926.                      }
  928.                      /* combine tests to reduce mispredicts -cfj */
  929.                      b1 = q_off_t <= H_query_etrim;
  930.                      b2 = s_off_t <= H_sub_etrim;
  931.                      
  932.                      if(0) if(H2_helper->maxsum1<=H_hsp_sum)break;
  933.                      
  934.                      if (!(b0|b1|b2) )
  935.                      {
  936.                         H2 = H2_helper->ptr;
  937.                         
  938.                         H_hsp_num=H2->hsp_link.num[index];
  939.                         H_hsp_sum=H2->hsp_link.sum[index];
  940.                         H_hsp_xsum=H2->hsp_link.xsum[index];
  941.                         H_hsp_link=H2;
  942.                      }
  943.                   } /* end for H2_index... */
  944.                } /* end if(H->score>cuttof[]) */
  945.                { 
  946.                   Int4 score=H->hsp->score;
  947.                   double new_xsum = 
  948.                      H_hsp_xsum + (score*(kbp[H->hsp->context]->Lambda));
  949.                   Int4 new_sum = H_hsp_sum + (score - cutoff[index]);
  950.                   
  951.                   H->hsp_link.sum[index] = new_sum;
  952.                   H->hsp_link.num[index] = H_hsp_num+1;
  953.                   H->hsp_link.link[index] = H_hsp_link;
  954.                   lh_helper[H_index].sum[index] = new_sum;
  955.                   lh_helper[H_index].maxsum1 = MAX(lh_helper[H_index-1].maxsum1, new_sum);
  956.                   /* Update this entry's 'next_larger' field */
  957.                   {
  958.                      Int4 cur_sum=lh_helper[H_index].sum[1];
  959.                      Int4 prev = H_index-1;
  960.                      Int4 prev_sum = lh_helper[prev].sum[1];
  961.                      while((cur_sum>=prev_sum) && (prev>0)){
  962.                         prev=lh_helper[prev].next_larger;
  963.                         prev_sum = lh_helper[prev].sum[1];
  964.                      }
  965.                      lh_helper[H_index].next_larger = prev;
  966.                   }
  967.                   
  968.                   if (new_sum >= maxscore) 
  969.                   {
  970.                      maxscore=new_sum;
  971.                      best[index]=H;
  972.                   }
  973.                   H->hsp_link.xsum[index] = new_xsum;
  974.                   if(H_hsp_link)
  975.                      ((LinkHSPStruct*)H_hsp_link)->linked_to++;
  976.                }
  977.             }
  978.             path_changed=0;
  979.             first_pass=0;
  980.          }
  981.          /********************************/
  982.          if (!ignore_small_gaps)
  983.          {
  984.             /* Select the best ordering method.
  985.                First we add back in the value cutoff[index] * the number
  986.                of links, as this was subtracted out for purposes of the
  987.                comparison above. */
  988.             best[0]->hsp_link.sum[0] +=
  989.                (best[0]->hsp_link.num[0])*cutoff[0];
  991.             prob[0] = BLAST_SmallGapSumE(kbp[query_context],
  992.                          gap_size,
  993.                          best[0]->hsp_link.num[0], best[0]->hsp_link.xsum[0],
  994.                          query_length, subject_length,
  995.                          BLAST_GapDecayDivisor(gap_decay_rate,
  996.                                               best[0]->hsp_link.num[0]) );
  998.             /* Adjust the e-value because we are performing multiple tests */
  999.             if( best[0]->hsp_link.num[0] > 1 ) {
  1000.               if( gap_prob == 0 || (prob[0] /= gap_prob) > INT4_MAX ) {
  1001.                 prob[0] = INT4_MAX;
  1002.               }
  1003.             }
  1005.             prob[1] = BLAST_LargeGapSumE(kbp[query_context],
  1006.                          best[1]->hsp_link.num[1],
  1007.                          best[1]->hsp_link.xsum[1],
  1008.                          query_length, subject_length,
  1009.                          BLAST_GapDecayDivisor(gap_decay_rate,
  1010.                                               best[1]->hsp_link.num[1]));
  1012.             if( best[1]->hsp_link.num[1] > 1 ) {
  1013.               if( 1 - gap_prob == 0 || (prob[1] /= 1 - gap_prob) > INT4_MAX ) {
  1014.                 prob[1] = INT4_MAX;
  1015.               }
  1016.             }
  1017.             ordering_method =
  1018.                prob[0]<=prob[1] ? BLAST_SMALL_GAPS : BLAST_LARGE_GAPS;
  1019.          }
  1020.          else
  1021.          {
  1022.             /* We only consider the case of big gaps. */
  1023.             best[1]->hsp_link.sum[1] +=
  1024.                (best[1]->hsp_link.num[1])*cutoff[1];
  1026.             prob[1] = BLAST_LargeGapSumE(kbp[query_context],
  1027.                          best[1]->hsp_link.num[1],
  1028.                          best[1]->hsp_link.xsum[1],
  1029.                          query_length, subject_length,
  1030.                          BLAST_GapDecayDivisor(gap_decay_rate,
  1031.                                               best[1]->hsp_link.num[1]));
  1032.             ordering_method = BLAST_LARGE_GAPS;
  1033.          }
  1035.          best[ordering_method]->start_of_chain = TRUE;
  1036.          
  1037.          /* AM: Support for query concatenation. */
  1038.          prob[ordering_method] *= 
  1039.             ((double)query_info->eff_searchsp_array[query_context] /
  1040.              ((double)subject_length*query_length));
  1041.          
  1042.          best[ordering_method]->hsp->evalue = prob[ordering_method];
  1043.          
  1044.          /* remove the links that have been ordered already. */
  1045.          if (best[ordering_method]->hsp_link.link[ordering_method])
  1046.             linked_set = TRUE;
  1047.          else
  1048.             linked_set = FALSE;
  1050.          if (best[ordering_method]->linked_to>0) path_changed=1;
  1051.          for (H=best[ordering_method]; H!=NULL;
  1052.               H=H->hsp_link.link[ordering_method]) 
  1053.          {
  1054.             if (H->linked_to>1) path_changed=1;
  1055.             H->linked_to=-1000;
  1056.             H->hsp_link.changed=1;
  1057.             /* record whether this is part of a linked set. */
  1058.             H->linked_set = linked_set;
  1059.             H->ordering_method = ordering_method;
  1060.             H->hsp->evalue = prob[ordering_method];
  1061.             if (H->next)
  1062.                (H->next)->prev=H->prev;
  1063.             if (H->prev)
  1064.                (H->prev)->next=H->next;
  1065.             number_of_hsps--;
  1066.          }
  1067.          
  1068.       } /* end while num_hsps... */
  1069. } /* end for frame_index ... */
  1071.    sfree(hp_frame_start);
  1072.    sfree(hp_frame_number);
  1073.    sfree(hp_start->hsp);
  1074.    sfree(hp_start);
  1076.    if (translated_query) {
  1077.       qsort(link_hsp_array,total_number_of_hsps,sizeof(LinkHSPStruct*), 
  1078.             rev_compare_hsps_transl);
  1079.       qsort(link_hsp_array, total_number_of_hsps,sizeof(LinkHSPStruct*), 
  1080.             fwd_compare_hsps_transl);
  1081.    } else {
  1082.       qsort(link_hsp_array,total_number_of_hsps,sizeof(LinkHSPStruct*), 
  1083.             rev_compare_hsps);
  1084.       qsort(link_hsp_array, total_number_of_hsps,sizeof(LinkHSPStruct*), 
  1085.             fwd_compare_hsps);
  1086.    }
  1088.    /* Sort by starting position. */
  1089.    
  1091. for (index=0, last_hsp=NULL;index<total_number_of_hsps; index++) 
  1092. {
  1093. H = link_hsp_array[index];
  1094. H->prev = NULL;
  1095. H->next = NULL;
  1096. }
  1098.    /* hook up the HSP's. */
  1099. first_hsp = NULL;
  1100. for (index=0, last_hsp=NULL;index<total_number_of_hsps; index++) 
  1101.    {
  1102. H = link_hsp_array[index];
  1104.       /* If this is not a single piece or the start of a chain, then Skip it. */
  1105.       if (H->linked_set == TRUE && H->start_of_chain == FALSE)
  1106. continue;
  1108.       /* If the HSP has no "link" connect the "next", otherwise follow the "link"
  1109.          chain down, connecting them with "next" and "prev". */
  1110. if (last_hsp == NULL)
  1111. first_hsp = H;
  1112. H->prev = last_hsp;
  1113. ordering_method = H->ordering_method;
  1114. if (H->hsp_link.link[ordering_method] == NULL)
  1115. {
  1116.          /* Grab the next HSP that is not part of a chain or the start of a chain */
  1117.          /* The "next" pointers are not hooked up yet in HSP's further down array. */
  1118.          index1=index;
  1119.          H2 = index1<(total_number_of_hsps-1) ? link_hsp_array[index1+1] : NULL;
  1120.          while (H2 && H2->linked_set == TRUE && 
  1121.                 H2->start_of_chain == FALSE)
  1122.          {
  1123.             index1++;
  1124.       H2 = index1<(total_number_of_hsps-1) ? link_hsp_array[index1+1] : NULL;
  1125.          }
  1126.          H->next= H2;
  1127. }
  1128. else
  1129. {
  1130. /* The first one has the number of links correct. */
  1131. num_links = H->hsp_link.num[ordering_method];
  1132. link = H->hsp_link.link[ordering_method];
  1133. while (link)
  1134. {
  1135. H->hsp->num = num_links;
  1136. H->sumscore = H->hsp_link.sum[ordering_method];
  1137. H->next = (LinkHSPStruct*) link;
  1138. H->prev = last_hsp;
  1139. last_hsp = H;
  1140. H = H->next;
  1141. if (H != NULL)
  1142.     link = H->hsp_link.link[ordering_method];
  1143. else
  1144.     break;
  1145. }
  1146. /* Set these for last link in chain. */
  1147. H->hsp->num = num_links;
  1148. H->sumscore = H->hsp_link.sum[ordering_method];
  1149.          /* Grab the next HSP that is not part of a chain or the start of a chain */
  1150.          index1=index;
  1151.          H2 = index1<(total_number_of_hsps-1) ? link_hsp_array[index1+1] : NULL;
  1152.          while (H2 && H2->linked_set == TRUE && 
  1153.                 H2->start_of_chain == FALSE)
  1154.    {
  1155.             index1++;
  1156.             H2 = index1<(total_number_of_hsps-1) ? link_hsp_array[index1+1] : NULL;
  1157. }
  1158.          H->next= H2;
  1159. H->prev = last_hsp;
  1160. }
  1161. last_hsp = H;
  1162. }
  1164.    /* The HSP's may be in a different order than they were before, 
  1165.       but first_hsp contains the first one. */
  1166.    for (index = 0, H = first_hsp; index < hsp_list->hspcnt; index++) {
  1167.       hsp_list->hsp_array[index] = H->hsp;
  1168.       /* Free the wrapper structure */
  1169.       H2 = H->next;
  1170.       sfree(H);
  1171.       H = H2;
  1172.    }
  1173.    sfree(link_hsp_array);
  1174.    sfree(lh_helper);
  1176.    return 0;
  1177. }
  1179. static void ConnectLinkHSPStructs(LinkHSPStruct** linkhsp_array, Int4 hspcnt,
  1180.                                   Uint1* subject_seq)
  1181. {
  1182.    Int4 index, index1, i;
  1183.    LinkHSPStruct* hsp;
  1185.    qsort(linkhsp_array, hspcnt, sizeof(LinkHSPStruct*), 
  1186.          sumscore_compare_hsps);
  1188.    hsp = linkhsp_array[0];
  1189.    for (index=0; index<hspcnt; hsp = hsp->next) {
  1190.       if (hsp->linked_set) {
  1191.          index1 = hsp->hsp->num;
  1192.          for (i=1; i < index1; i++, hsp = hsp->next) {
  1193.             hsp->next->hsp->evalue = hsp->hsp->evalue; 
  1194.             hsp->next->hsp->num = hsp->hsp->num;
  1195.             hsp->next->sumscore = hsp->sumscore;
  1196.             if (FindSpliceJunction(subject_seq, hsp->hsp, hsp->next->hsp)) {
  1197.                /* Kludge: ordering_method here would indicate existence of
  1198.                   splice junctions(s) */
  1199.                hsp->hsp->splice_junction++;
  1200.                hsp->next->hsp->splice_junction++;
  1201.             } else {
  1202.                hsp->hsp->splice_junction--;
  1203.                hsp->next->hsp->splice_junction--;
  1204.             }
  1205.          }
  1206.       } 
  1207.       while (++index < hspcnt)
  1208.          if (!linkhsp_array[index]->linked_set ||
  1209.              linkhsp_array[index]->start_of_chain)
  1210.             break;
  1211.       if (index == hspcnt) {
  1212.          hsp->next = NULL;
  1213.          break;
  1214.       }
  1215.       hsp->next = linkhsp_array[index];
  1216.    }
  1217. }
  1219. static Int2
  1220. new_link_hsps(Uint1 program_number, BlastHSPList* hsp_list, 
  1221.    BlastQueryInfo* query_info, BLAST_SequenceBlk* subject,
  1222.    BlastScoreBlk* sbp, const BlastHitSavingParameters* hit_params)
  1223. {
  1224.    BlastHSP** hsp_array;
  1225.    LinkHSPStruct** score_hsp_array,** offset_hsp_array,** end_hsp_array;
  1226.    LinkHSPStruct* hsp,* head_hsp,* best_hsp,* var_hsp;
  1227.    Int4 hspcnt, index, index1, i;
  1228.    double best_evalue, evalue;
  1229.    Int4 sumscore, best_sumscore = 0;
  1230.    Boolean reverse_link;
  1231.    Int4 longest_intron = hit_params->options->longest_intron;
  1232.    LinkHSPStruct** link_hsp_array;
  1234.    hspcnt = hsp_list->hspcnt;
  1235.    hsp_array = hsp_list->hsp_array;
  1237.    link_hsp_array = (LinkHSPStruct**) malloc(hspcnt*sizeof(LinkHSPStruct*));
  1238.    for (index = 0; index < hspcnt; ++index) {
  1239.       link_hsp_array[index] = (LinkHSPStruct*) calloc(1, sizeof(LinkHSPStruct));
  1240.       link_hsp_array[index]->hsp = hsp_array[index];
  1241.       link_hsp_array[index]->linked_set = FALSE;
  1242.       hsp_array[index]->num = 1;
  1243.       hsp_array[index]->splice_junction = FALSE;
  1244.    }
  1246.    score_hsp_array = (LinkHSPStruct**) malloc(hspcnt*sizeof(LinkHSPStruct*));
  1247.    offset_hsp_array = (LinkHSPStruct**) malloc(hspcnt*sizeof(LinkHSPStruct*));
  1248.    end_hsp_array = (LinkHSPStruct**) malloc(hspcnt*sizeof(LinkHSPStruct*));
  1250.    memcpy(score_hsp_array, link_hsp_array, hspcnt*sizeof(LinkHSPStruct*));
  1251.    memcpy(offset_hsp_array, link_hsp_array, hspcnt*sizeof(LinkHSPStruct*));
  1252.    memcpy(end_hsp_array, link_hsp_array, hspcnt*sizeof(LinkHSPStruct*));
  1253.    qsort(offset_hsp_array, hspcnt, sizeof(LinkHSPStruct*), fwd_compare_hsps);
  1254.    qsort(end_hsp_array, hspcnt, sizeof(LinkHSPStruct*), end_compare_hsps);
  1256.    qsort(score_hsp_array, hspcnt, sizeof(LinkHSPStruct*), 
  1257.             sumscore_compare_hsps);
  1258.       
  1259.    head_hsp = NULL;
  1260.    for (index = 0; index < hspcnt && score_hsp_array[index]; ) {
  1261.       if (!head_hsp) {
  1262.          while (index<hspcnt && score_hsp_array[index] && 
  1263.                 score_hsp_array[index]->linked_set)
  1264.             index++;
  1265.          if (index==hspcnt)
  1266.             break;
  1267.          head_hsp = score_hsp_array[index];
  1268.       }
  1269.       best_evalue = head_hsp->hsp->evalue;
  1270.       best_hsp = NULL;
  1271.       reverse_link = FALSE;
  1272.       
  1273.       if (head_hsp->linked_set)
  1274.          for (var_hsp = head_hsp, i=1; i<head_hsp->hsp->num; 
  1275.               var_hsp = var_hsp->next, i++);
  1276.       else
  1277.          var_hsp = head_hsp;
  1279.       index1 = hsp_binary_search(offset_hsp_array, hspcnt,
  1280.                                  var_hsp->hsp->query.end - WINDOW_SIZE, TRUE);
  1281.       while (index1 < hspcnt && offset_hsp_array[index1]->hsp->query.offset < 
  1282.              var_hsp->hsp->query.end + WINDOW_SIZE) {
  1283.          hsp = offset_hsp_array[index1++];
  1284.          /* If this is already part of a linked set, disregard it */
  1285.          if (hsp == var_hsp || hsp == head_hsp || 
  1286.              (hsp->linked_set && !hsp->start_of_chain))
  1287.             continue;
  1288.          /* Check if the subject coordinates are consistent with query */
  1289.          if ((hsp->hsp->subject.offset < 
  1290.              var_hsp->hsp->subject.end - WINDOW_SIZE) || 
  1291.              (hsp->hsp->subject.offset > 
  1292.               var_hsp->hsp->subject.end + longest_intron))
  1293.             continue;
  1294.          /* Check if the e-value for the combined two HSPs is better than for
  1295.             one of them */
  1296.          if ((evalue = SumHSPEvalue(program_number, sbp, query_info, subject, 
  1297.                                     hit_params, head_hsp, hsp, &sumscore)) < 
  1298.              MIN(best_evalue, hsp->hsp->evalue)) {
  1299.             best_hsp = hsp;
  1300.             best_evalue = evalue;
  1301.             best_sumscore = sumscore;
  1302.          }
  1303.       }
  1304.       index1 = hsp_binary_search(end_hsp_array, hspcnt,
  1305.                                  head_hsp->hsp->query.offset - WINDOW_SIZE, FALSE);
  1306.       while (index1 < hspcnt && end_hsp_array[index1]->hsp->query.end < 
  1307.              head_hsp->hsp->query.offset + WINDOW_SIZE) {
  1308.          hsp = end_hsp_array[index1++];
  1310.          /* Check if the subject coordinates are consistent with query */
  1311.          if (hsp == head_hsp || 
  1312.              hsp->hsp->subject.end > head_hsp->hsp->subject.offset + WINDOW_SIZE || 
  1313.              hsp->hsp->subject.end < head_hsp->hsp->subject.offset - longest_intron)
  1314.             continue;
  1315.          if (hsp->linked_set) {
  1316.             for (var_hsp = hsp, i=1; var_hsp->start_of_chain == FALSE; 
  1317.                  var_hsp = var_hsp->prev, i++);
  1318.             if (i<var_hsp->hsp->num || var_hsp == head_hsp)
  1319.                continue;
  1320.          } else
  1321.             var_hsp = hsp;
  1322.          /* Check if the e-value for the combined two HSPs is better than for
  1323.             one of them */
  1324.          if ((evalue = SumHSPEvalue(program_number, sbp, query_info, subject, 
  1325.                           hit_params, var_hsp, head_hsp, &sumscore)) < 
  1326.              MIN(var_hsp->hsp->evalue, best_evalue)) {
  1327.             best_hsp = hsp;
  1328.             best_evalue = evalue;
  1329.             best_sumscore = sumscore;
  1330.             reverse_link = TRUE;
  1331.          }
  1332.       }
  1333.          
  1334.       /* Link these HSPs together */
  1335.       if (best_hsp != NULL) {
  1336.          if (!reverse_link) {
  1337.             head_hsp->start_of_chain = TRUE;
  1338.             head_hsp->sumscore = best_sumscore;
  1339.             head_hsp->hsp->evalue = best_evalue;
  1340.             best_hsp->start_of_chain = FALSE;
  1341.             if (head_hsp->linked_set) 
  1342.                for (var_hsp = head_hsp, i=1; i<head_hsp->hsp->num; 
  1343.                     var_hsp = var_hsp->next, i++);
  1344.             else 
  1345.                var_hsp = head_hsp;
  1346.             var_hsp->next = best_hsp;
  1347.             best_hsp->prev = var_hsp;
  1348.             head_hsp->hsp->num += best_hsp->hsp->num;
  1349.          } else {
  1350.             best_hsp->next = head_hsp;
  1351.             head_hsp->prev = best_hsp;
  1352.             if (best_hsp->linked_set) {
  1353.                for (var_hsp = best_hsp; 
  1354.                     var_hsp->start_of_chain == FALSE; 
  1355.                     var_hsp = var_hsp->prev);
  1356.             } else
  1357.                   var_hsp = best_hsp;
  1358.             var_hsp->start_of_chain = TRUE;
  1359.             var_hsp->sumscore = best_sumscore;
  1360.             var_hsp->hsp->evalue = best_evalue;
  1361.             var_hsp->hsp->num += head_hsp->hsp->num;
  1362.             head_hsp->start_of_chain = FALSE;
  1363.          }
  1364.          
  1365.          head_hsp->linked_set = best_hsp->linked_set = TRUE;
  1366.          if (reverse_link)
  1367.             head_hsp = var_hsp;
  1368.       } else {
  1369.          head_hsp = NULL;
  1370.          index++;
  1371.       }
  1372.    }
  1373.   
  1374.    ConnectLinkHSPStructs(score_hsp_array, hspcnt, 
  1375.                          subject->sequence_start);
  1376.    head_hsp = score_hsp_array[0];
  1378.    sfree(score_hsp_array);
  1379.    sfree(offset_hsp_array);
  1380.    sfree(end_hsp_array);
  1382.    /* Place HSPs in the original HSP array in their new order */
  1383.    for (index=0; head_hsp && index<hspcnt; index++) {
  1384.       hsp_list->hsp_array[index] = head_hsp->hsp;
  1385.       var_hsp = head_hsp->next;
  1386.       sfree(head_hsp);
  1387.       head_hsp = var_hsp;
  1388.    }
  1390.    sfree(link_hsp_array);
  1391.    return 0;
  1392. }
  1394. Int2 
  1395. BLAST_LinkHsps(Uint1 program_number, BlastHSPList* hsp_list, 
  1396.    BlastQueryInfo* query_info, BLAST_SequenceBlk* subject, 
  1397.    BlastScoreBlk* sbp, const BlastHitSavingParameters* hit_params,
  1398.    Boolean gapped_calculation)
  1399. {
  1400. if (hsp_list && hsp_list->hspcnt > 0)
  1401. {
  1402.       /* Link up the HSP's for this hsp_list. */
  1403.       if (program_number != blast_type_tblastn ||
  1404.           hit_params->options->longest_intron <= 0)
  1405.       {
  1406.          link_hsps(program_number, hsp_list, query_info, subject, sbp, 
  1407.                    hit_params, gapped_calculation);
  1408.          /* The HSP's may be in a different order than they were before, 
  1409.             but hsp contains the first one. */
  1410.       } else {
  1411.          /* Calculate individual HSP e-values first - they'll be needed to
  1412.             compare with sum e-values. */
  1413.          Blast_HSPListGetEvalues(program_number, query_info, 
  1414.                                     hsp_list, gapped_calculation, sbp);
  1415.          
  1416.          new_link_hsps(program_number, hsp_list, query_info, subject, sbp, 
  1417.                        hit_params);
  1418.       }
  1419. }
  1421. return 0;
  1422. }