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

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

生物技术

开发平台：

C/C++

blast_filter.c：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: blast_filter.c,v $
  4.  * PRODUCTION Revision 1000.3  2004/06/01 18:07:16  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.45
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /* $Id: blast_filter.c,v 1000.3 2004/06/01 18:07:16 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 blast_filter.c
  40.  * All code related to query sequence masking/filtering for BLAST
  41.  */
  43. static char const rcsid[] = 
  44.     "$Id: blast_filter.c,v 1000.3 2004/06/01 18:07:16 gouriano Exp $";
  46. #include <algo/blast/core/blast_def.h>
  47. #include <algo/blast/core/blast_util.h>
  48. #include <algo/blast/core/blast_filter.h>
  49. #include <algo/blast/core/blast_dust.h>
  50. #include <algo/blast/core/blast_seg.h>
  51. #ifdef CC_FILTER_ALLOWED
  52. #include <algo/blast/core/urkpcc.h>
  53. #endif
  55. /* The following function will replace BlastSetUp_CreateSSeqRange */
  57. BlastSeqLoc* BlastSeqLocNew(Int4 from, Int4 to)
  58. {
  59.    BlastSeqLoc* loc = (BlastSeqLoc*) calloc(1, sizeof(BlastSeqLoc));
  60.    SSeqRange* di = (SSeqRange*) malloc(sizeof(SSeqRange));
  62.    di->left = from;
  63.    di->right = to;
  64.    loc->ptr = di;
  65.    return loc;
  66. }
  68. BlastSeqLoc* BlastSeqLocFree(BlastSeqLoc* loc)
  69. {
  70.    SSeqRange* dintp;
  71.    BlastSeqLoc* next_loc;
  73.    while (loc) {
  74.       next_loc = loc->next;
  75.       dintp = (SSeqRange*) loc->ptr;
  76.       sfree(dintp);
  77.       sfree(loc);
  78.       loc = next_loc;
  79.    }
  80.    return NULL;
  81. }
  83. BlastMaskLoc* BlastMaskLocNew(Int4 index, BlastSeqLoc *loc_list)
  84. {
  85.       BlastMaskLoc* retval = (BlastMaskLoc *) calloc(1, sizeof(BlastMaskLoc));
  86.       retval->index = index;
  87.       retval->loc_list = loc_list;
  88.       return retval;
  89. }
  91. BlastMaskLoc* BlastMaskLocFree(BlastMaskLoc* mask_loc)
  92. {
  93.    BlastMaskLoc* next_loc;
  94.    while (mask_loc) {
  95.       next_loc = mask_loc->next;
  96.       BlastSeqLocFree(mask_loc->loc_list);
  97.       sfree(mask_loc);
  98.       mask_loc = next_loc;
  99.    }
  100.    return NULL;
  101. }
  103. /** Used for qsort, compares two SeqLoc's by starting position. */
  104. static int SSeqRangeSortByStartPosition(const void *vp1, const void *vp2)
  105. {
  106.    ListNode* v1 = *((ListNode**) vp1);
  107.    ListNode* v2 = *((ListNode**) vp2);
  108.    SSeqRange* loc1 = (SSeqRange*) v1->ptr;
  109.    SSeqRange* loc2 = (SSeqRange*) v2->ptr;
  110.    
  111.    if (loc1->left < loc2->left)
  112.       return -1;
  113.    else if (loc1->left > loc2->left)
  114.       return 1;
  115.    else
  116.       return 0;
  117. }
  119. /* This will go in place of CombineSeqLocs to combine filtered locations */
  120. Int2
  121. CombineMaskLocations(BlastSeqLoc* mask_loc, BlastSeqLoc* *mask_loc_out,
  122.                      Int4 link_value)
  123. {
  124.    Int2 status=0; /* return value. */
  125.    Int4 start, stop; /* USed to merge overlapping SeqLoc's. */
  126.    SSeqRange* di = NULL,* di_next = NULL,* di_tmp=NULL;
  127.    BlastSeqLoc* loc_head=NULL,* last_loc=NULL,* loc_var=NULL;
  128.    BlastSeqLoc* new_loc = NULL,* new_loc_last = NULL;
  129.    
  130.    if (!mask_loc) {
  131.       *mask_loc_out = NULL;
  132.       return 0;
  133.    }
  135.    /* Put all the SeqLoc's into one big linked list. */
  136.    loc_head = last_loc = 
  137.       (BlastSeqLoc*) BlastMemDup(mask_loc, sizeof(BlastSeqLoc));
  138.          
  139.    /* Copy all locations, so loc points at the end of the chain */
  140.    while (last_loc->next) {
  141.       last_loc->next = 
  142.          (BlastSeqLoc*) BlastMemDup(last_loc->next, sizeof(BlastSeqLoc));
  143.       last_loc = last_loc->next;
  144.    }
  145.    
  146.    /* Sort them by starting position. */
  147.    loc_head = (BlastSeqLoc*) 
  148.       ListNodeSort ((ListNode*) loc_head, 
  149.                    SSeqRangeSortByStartPosition);
  150.    
  151.    di = (SSeqRange*) loc_head->ptr;
  152.    start = di->left;
  153.    stop = di->right;
  154.    loc_var = loc_head;
  155.    
  156.    while (loc_var) {
  157.       di = loc_var->ptr;
  158.       if (loc_var->next)
  159.          di_next = loc_var->next->ptr;
  160.       if (di_next && ((stop + link_value) > di_next->left)) {
  161.          stop = MAX(stop, di_next->right);
  162.       } else {
  163.          di_tmp = (SSeqRange*) malloc(sizeof(SSeqRange));
  164.          di_tmp->left = start;
  165.          di_tmp->right = stop;
  166.          if (!new_loc)
  167.             new_loc_last = ListNodeAddPointer(&new_loc, 0, di_tmp);
  168.          else
  169.             new_loc_last = ListNodeAddPointer(&new_loc_last, 0, di_tmp);
  170.          if (loc_var->next) {
  171.              start = di_next->left;
  172.              stop = di_next->right;
  173.          }
  174.       }
  175.       loc_var = loc_var->next;
  176.       di_next = NULL;
  177.    }
  179.    *mask_loc_out = new_loc;
  180.       
  181.    /* Free memory allocated for the temporary list of SeqLocs */
  182.    while (loc_head) {
  183.       loc_var = loc_head->next;
  184.       sfree(loc_head);
  185.       loc_head = loc_var;
  186.    }
  187.    return status;
  188. }
  190. Int2 
  191. BLAST_ComplementMaskLocations(Uint1 program_number, 
  192.    BlastQueryInfo* query_info, 
  193.    BlastMaskLoc* mask_loc, BlastSeqLoc* *complement_mask) 
  194. {
  195.    Int4 start_offset, end_offset, filter_start, filter_end;
  196.    Int4 context, index;
  197.    BlastSeqLoc* loc,* last_loc = NULL,* start_loc = NULL;
  198.    SSeqRange* double_int = NULL,* di;
  199.    Boolean first; /* Specifies beginning of query. */
  200.    Boolean last_interval_open=TRUE; /* if TRUE last interval needs to be closed. */
  201.    Boolean reverse = FALSE;
  202.    const Boolean k_is_na = (program_number == blast_type_blastn);
  204.    if (complement_mask == NULL)
  205. return -1;
  207.    *complement_mask = NULL;
  209.    for (context = query_info->first_context; 
  210.         context <= query_info->last_context; ++context) {
  211.       start_offset = query_info->context_offsets[context];
  212.       end_offset = query_info->context_offsets[context+1] - 2;
  213.       /* For blastn: check if this strand is not searched at all */
  214.       if (end_offset < start_offset)
  215.           continue;
  216.       index = (k_is_na ? context / 2 : context);
  217.       reverse = (k_is_na && ((context & 1) != 0));
  218.       first = TRUE;
  220.       if (!reverse) {
  221.          for ( ; mask_loc && mask_loc->index < index; 
  222.                mask_loc = mask_loc->next);
  223.       }
  225.       if (!mask_loc || (mask_loc->index > index) ||
  226.           !mask_loc->loc_list) {
  227.          /* No masks for this context */
  228.          double_int = (SSeqRange*) calloc(1, sizeof(SSeqRange));
  229.          double_int->left = start_offset;
  230.          double_int->right = end_offset;
  231.          if (!last_loc)
  232.             last_loc = ListNodeAddPointer(complement_mask, 0, double_int);
  233.          else 
  234.             last_loc = ListNodeAddPointer(&last_loc, 0, double_int);
  235.          continue;
  236.       }
  237.       
  238.       if (reverse) {
  239.          BlastSeqLoc* prev_loc = NULL;
  240.          /* Reverse the order of the locations */
  241.          for (start_loc = mask_loc->loc_list; start_loc; 
  242.               start_loc = start_loc->next) {
  243.             loc = (BlastSeqLoc*) BlastMemDup(start_loc, sizeof(BlastSeqLoc));
  244.             loc->next = prev_loc;
  245.             prev_loc = loc;
  246.          }
  247.          /* Save where this list starts, so it can be freed later */
  248.          start_loc = loc;
  249.       } else {
  250.          loc = mask_loc->loc_list;
  251.       }
  253.       for ( ; loc; loc = loc->next) {
  254.          di = loc->ptr;
  255.          if (reverse) {
  256.             filter_start = end_offset - di->right;
  257.             filter_end = end_offset - di->left;
  258.          } else {
  259.             filter_start = start_offset + di->left;
  260.             filter_end = start_offset + di->right;
  261.          }
  262.          /* The canonical "state" at the top of this 
  263.             while loop is that a SSeqRange has been 
  264.             created and one field was filled in on the 
  265.             last iteration. The first time this loop is 
  266.             entered in a call to the funciton this is not
  267.             true and the following "if" statement moves 
  268.             everything to the canonical state. */
  269.          if (first) {
  270.             last_interval_open = TRUE;
  271.             first = FALSE;
  272.             double_int = (SSeqRange*) calloc(1, sizeof(SSeqRange));
  273.             
  274.             if (filter_start > start_offset) {
  275.                /* beginning of sequence not filtered */
  276.                double_int->left = start_offset;
  277.             } else {
  278.                /* beginning of sequence filtered */
  279.                double_int->left = filter_end + 1;
  280.                continue;
  281.             }
  282.          }
  284.          double_int->right = filter_start - 1;
  286.          if (!last_loc)
  287.             last_loc = ListNodeAddPointer(complement_mask, 0, double_int);
  288.          else 
  289.             last_loc = ListNodeAddPointer(&last_loc, 0, double_int);
  290.          if (filter_end >= end_offset) {
  291.             /* last masked region at end of sequence */
  292.             last_interval_open = FALSE;
  293.             break;
  294.          } else {
  295.             double_int = (SSeqRange*) calloc(1, sizeof(SSeqRange));
  296.                double_int->left = filter_end + 1;
  297.          }
  298.       }
  300.       if (reverse) {
  301.          start_loc = ListNodeFree(start_loc);
  302.       }
  303.       
  304.       if (last_interval_open) {
  305.          /* Need to finish SSeqRange* for last interval. */
  306.          double_int->right = end_offset;
  307.          if (!last_loc)
  308.             last_loc = ListNodeAddPointer(complement_mask, 0, double_int);
  309.          else 
  310.             last_loc = ListNodeAddPointer(&last_loc, 0, double_int);
  311.          last_interval_open = FALSE;
  312.       }
  313.    }
  314.    return 0;
  315. }
  317. #define BLASTOPTIONS_BUFFER_SIZE 128
  319. /** Parses options used for dust.
  320.  * @param ptr buffer containing instructions. [in]
  321.  * @param level sets level for dust. [out]
  322.  * @param window sets window for dust [out]
  323.  * @param cutoff sets cutoff for dust. [out] 
  324.  * @param linker sets linker for dust. [out]
  325. */
  326. static Int2
  327. parse_dust_options(const char *ptr, Int2* level, Int2* window,
  328. Int2* cutoff, Int2* linker)
  330. {
  331. char buffer[BLASTOPTIONS_BUFFER_SIZE];
  332. Int4 arg, index, index1, window_pri=-1, linker_pri=-1, level_pri=-1, cutoff_pri=-1;
  333. long tmplong;
  335. arg = 0;
  336. index1 = 0;
  337. for (index=0; index<BLASTOPTIONS_BUFFER_SIZE; index++)
  338. {
  339. if (*ptr == ' ' || *ptr == NULLB)
  340. {
  341. buffer[index1] = NULLB;
  342. index1 = 0;
  343. switch(arg) {
  344. case 0:
  345. sscanf(buffer, "%ld", &tmplong);
  346. level_pri = tmplong;
  347. break;
  348. case 1:
  349. sscanf(buffer, "%ld", &tmplong);
  350. window_pri = tmplong;
  351. break;
  352. case 2:
  353. sscanf(buffer, "%ld", &tmplong);
  354. cutoff_pri = tmplong;
  355. break;
  356. case 3:
  357. sscanf(buffer, "%ld", &tmplong);
  358. linker_pri = tmplong;
  359. break;
  360. default:
  361. break;
  362. }
  364. arg++;
  365. while (*ptr == ' ')
  366. ptr++;
  368. /* end of the buffer. */
  369. if (*ptr == NULLB)
  370. break;
  371. }
  372. else
  373. {
  374. buffer[index1] = *ptr; ptr++;
  375. index1++;
  376. }
  377. }
  379. *level = level_pri; 
  380. *window = window_pri; 
  381. *cutoff = cutoff_pri; 
  382. *linker = linker_pri; 
  384. return 0;
  385. }
  387. /** parses a string to set three seg options. 
  388.  * @param ptr buffer containing instructions [in]
  389.  * @param window returns "window" for seg algorithm. [out]
  390.  * @param locut returns "locut" for seg. [out]
  391.  * @param hicut returns "hicut" for seg. [out]
  392. */
  393. static Int2
  394. parse_seg_options(const char *ptr, Int4* window, double* locut, double* hicut)
  396. {
  397. char buffer[BLASTOPTIONS_BUFFER_SIZE];
  398. Int4 arg, index, index1; 
  399. long tmplong;
  400. double tmpdouble;
  402. arg = 0;
  403. index1 = 0;
  404. for (index=0; index<BLASTOPTIONS_BUFFER_SIZE; index++)
  405. {
  406. if (*ptr == ' ' || *ptr == NULLB)
  407. {
  408. buffer[index1] = NULLB;
  409. index1 = 0;
  410. switch(arg) {
  411. case 0:
  412. sscanf(buffer, "%ld", &tmplong);
  413. *window = tmplong;
  414. break;
  415. case 1:
  416. sscanf(buffer, "%le", &tmpdouble);
  417. *locut = tmpdouble;
  418. break;
  419. case 2:
  420. sscanf(buffer, "%le", &tmpdouble);
  421. *hicut = tmpdouble;
  422. break;
  423. default:
  424. break;
  425. }
  427. arg++;
  428. while (*ptr == ' ')
  429. ptr++;
  431. /* end of the buffer. */
  432. if (*ptr == NULLB)
  433. break;
  434. }
  435. else
  436. {
  437. buffer[index1] = *ptr; ptr++;
  438. index1++;
  439. }
  440. }
  442. return 0;
  443. }
  445. #ifdef CC_FILTER_ALLOWED
  446. /** Coiled-coiled algorithm parameters. 
  447.  * @param ptr buffer containing instructions. [in]
  448.  * @param window returns window for coil-coiled algorithm. [out]
  449.  * @param cutoff cutoff for coil-coiled algorithm [out]
  450.  * @param linker returns linker [out]
  451. */
  452. static Int2
  453. parse_cc_options(const char *ptr, Int4* window, double* cutoff, Int4* linker)
  455. {
  456. char buffer[BLASTOPTIONS_BUFFER_SIZE];
  457. Int4 arg, index, index1;
  458. long tmplong;
  459. double tmpdouble;
  461. arg = 0;
  462. index1 = 0;
  463. for (index=0; index<BLASTOPTIONS_BUFFER_SIZE; index++)
  464. {
  465. if (*ptr == ' ' || *ptr == NULLB)
  466. {
  467. buffer[index1] = NULLB;
  468. index1 = 0;
  469. switch(arg) {
  470. case 0:
  471. sscanf(buffer, "%ld", &tmplong);
  472. *window = tmplong;
  473. break;
  474. case 1:
  475. sscanf(buffer, "%le", &tmpdouble);
  476. *cutoff = tmpdouble;
  477. break;
  478. case 2:
  479. sscanf(buffer, "%ld", &tmplong);
  480. *linker = tmplong;
  481. break;
  482. default:
  483. break;
  484. }
  486. arg++;
  487. while (*ptr == ' ')
  488. ptr++;
  490. /* end of the buffer. */
  491. if (*ptr == NULLB)
  492. break;
  493. }
  494. else
  495. {
  496. buffer[index1] = *ptr; ptr++;
  497. index1++;
  498. }
  499. }
  501. return 0;
  502. }
  503. #endif
  505. /** Copies filtering commands for one filtering algorithm from "instructions" to
  506.  * "buffer". 
  507.  * ";" is a delimiter for the commands for different algorithms, so it stops
  508.  * copying when a ";" is found. 
  509.  * Example filtering string: "m L; R -d rodents.lib"
  510.  * @param instructions filtering commands [in] 
  511.  * @param buffer filled with filtering commands for one algorithm. [out]
  512. */
  513. static const char *
  514. BlastSetUp_load_options_to_buffer(const char *instructions, char* buffer)
  515. {
  516. Boolean not_started=TRUE;
  517. char* buffer_ptr;
  518. const char *ptr;
  519. Int4 index;
  521. ptr = instructions;
  522. buffer_ptr = buffer;
  523. for (index=0; index<BLASTOPTIONS_BUFFER_SIZE && *ptr != NULLB; index++)
  524. {
  525. if (*ptr == ';')
  526. { /* ";" is a delimiter for different filtering algorithms. */
  527. ptr++;
  528. break;
  529. }
  530. /* Remove blanks at the beginning. */
  531. if (not_started && *ptr == ' ')
  532. {
  533. ptr++;
  534. }
  535. else
  536. {
  537. not_started = FALSE;
  538. *buffer_ptr = *ptr;
  539. buffer_ptr++; ptr++;
  540. }
  541. }
  543. *buffer_ptr = NULLB;
  545. if (not_started == FALSE)
  546. { /* Remove trailing blanks. */
  547. buffer_ptr--;
  548. while (*buffer_ptr == ' ' && buffer_ptr > buffer)
  549. {
  550. *buffer_ptr = NULLB;
  551. buffer_ptr--;
  552. }
  553. }
  555. return ptr;
  556. }
  558. #define CC_WINDOW 22
  559. #define CC_CUTOFF 40.0
  560. #define CC_LINKER 32
  562. Int2
  563. BlastSetUp_Filter(Uint1 program_number, Uint1* sequence, Int4 length, 
  564.    Int4 offset, const char* instructions, Boolean *mask_at_hash, 
  565.    BlastSeqLoc* *seqloc_retval)
  566. {
  567. Boolean do_default=FALSE, do_seg=FALSE, do_dust=FALSE; 
  568. char* buffer=NULL;
  569.    const char *ptr;
  570. Int2 seqloc_num;
  571. Int2 status=0; /* return value. */
  572. Int2 window_dust, level_dust, minwin_dust, linker_dust;
  573.    BlastSeqLoc* dust_loc = NULL,* seg_loc = NULL;
  574. SegParameters* sparamsp=NULL;
  575. #ifdef CC_FILTER_ALLOWED
  576.    Boolean do_coil_coil = FALSE;
  577.    BlastSeqLoc* cc_loc = NULL;
  578. PccDatPtr pccp;
  579.    Int4 window_cc, linker_cc;
  580. double cutoff_cc;
  581. #endif
  582. #ifdef TEMP_BLAST_OPTIONS
  583.    /* TEMP_BLAST_OPTIONS is set to zero until these are implemented. */
  584.    BlastSeqLoc* vs_loc = NULL,* repeat_loc = NULL;
  585. BLAST_OptionsBlkPtr repeat_options, vs_options;
  586. Boolean do_repeats=FALSE;  /* screen for orgn. specific repeats. */
  587. Boolean do_vecscreen=FALSE; /* screen for vector contamination. */
  588. Boolean myslp_allocated;
  589. char* repeat_database=NULL,* vs_database=NULL,* error_msg;
  590. SeqLocPtr repeat_slp=NULL, vs_slp=NULL;
  591. SeqAlignPtr seqalign;
  592. SeqLocPtr myslp, seqloc_var, seqloc_tmp;
  593. ListNode* vnp=NULL,* vnp_var;
  594. #endif
  596. #ifdef CC_FILTER_ALLOWED
  597. cutoff_cc = CC_CUTOFF;
  598. #endif
  600.    if (!seqloc_retval) 
  601.       return -1;
  603. /* FALSE is the default right now. */
  604. if (mask_at_hash)
  605.       *mask_at_hash = FALSE;
  606.    
  607.    *seqloc_retval = NULL;
  609. if (instructions == NULL || strcasecmp(instructions, "F") == 0)
  610. return status;
  612. /* parameters for dust. */
  613. /* -1 indicates defaults. */
  614. level_dust = -1;
  615. window_dust = -1;
  616. minwin_dust = -1;
  617. linker_dust = -1;
  618. if (strcasecmp(instructions, "T") == 0)
  619. { /* do_default actually means seg for proteins and dust for nt. */
  620. do_default = TRUE;
  621. }
  622. else
  623. {
  624. buffer = (char*) calloc(strlen(instructions), sizeof(char));
  625. ptr = instructions;
  626. /* allow old-style filters when m cannot be followed by the ';' */
  627. if (*ptr == 'm' && ptr[1] == ' ')
  628. {
  629. if (mask_at_hash)
  630. *mask_at_hash = TRUE;
  631. ptr += 2;
  632. }
  633. while (*ptr != NULLB)
  634. {
  635. if (*ptr == 'S')
  636. {
  637. sparamsp = SegParametersNewAa();
  638. sparamsp->overlaps = TRUE; /* merge overlapping segments. */
  639. ptr = BlastSetUp_load_options_to_buffer(ptr+1, buffer);
  640. if (buffer[0] != NULLB)
  641. {
  642. parse_seg_options(buffer, &sparamsp->window, &sparamsp->locut, &sparamsp->hicut);
  643. }
  644. do_seg = TRUE;
  645. }
  646. else if (*ptr == 'C')
  647. {
  648. #ifdef CC_FILTER_ALLOWED
  649. ptr = BlastSetUp_load_options_to_buffer(ptr+1, buffer);
  650. window_cc = CC_WINDOW;
  651. cutoff_cc = CC_CUTOFF;
  652. linker_cc = CC_LINKER;
  653. if (buffer[0] != NULLB)
  654. parse_cc_options(buffer, &window_cc, &cutoff_cc, &linker_cc);
  655. do_coil_coil = TRUE;
  656. #endif
  657. }
  658. else if (*ptr == 'D')
  659. {
  660. ptr = BlastSetUp_load_options_to_buffer(ptr+1, buffer);
  661. if (buffer[0] != NULLB)
  662. parse_dust_options(buffer, &level_dust, &window_dust, &minwin_dust, &linker_dust);
  663. do_dust = TRUE;
  664. }
  665. #ifdef TEMP_BLAST_OPTIONS
  666. else if (*ptr == 'R')
  667. {
  668. repeat_options = BLASTOptionNew("blastn", TRUE);
  669. repeat_options->expect_value = 0.1;
  670. repeat_options->penalty = -1;
  671. repeat_options->wordsize = 11;
  672. repeat_options->gap_x_dropoff_final = 90;
  673. repeat_options->dropoff_2nd_pass = 40;
  674. repeat_options->gap_open = 2;
  675. repeat_options->gap_extend = 1;
  676. ptr = BlastSetUp_load_options_to_buffer(ptr+1, buffer);
  677. if (buffer[0] != NULLB)
  678.                                    parse_blast_options(repeat_options,
  679.                                       buffer, &error_msg, &repeat_database,
  680.                                       NULL, NULL);
  681. if (repeat_database == NULL)
  682.                                    repeat_database = strdup("humlines.lib humsines.lib retrovir.lib");
  683. do_repeats = TRUE;
  684. }
  685. else if (*ptr == 'V')
  686. {
  687. vs_options = VSBlastOptionNew();
  688. ptr = BlastSetUp_load_options_to_buffer(ptr+1, buffer);
  689. if (buffer[0] != NULLB)
  690.                                    parse_blast_options(vs_options, buffer,
  691.                                       &error_msg, &vs_database, NULL, NULL); 
  692. vs_options = BLASTOptionDelete(vs_options);
  693. if (vs_database == NULL)
  694.                                    vs_database = strdup("UniVec_Core");
  695. do_vecscreen = TRUE;
  696. }
  697. #endif
  698. else if (*ptr == 'L')
  699. { /* do low-complexity filtering; dust for blastn, otherwise seg.*/
  700. do_default = TRUE;
  701. ptr++;
  702. }
  703. else if (*ptr == 'm')
  704. {
  705. if (mask_at_hash)
  706. *mask_at_hash = TRUE;
  707. ptr++;
  708. }
  709. else
  710. { /* Nothing applied. */
  711. ptr++;
  712. }
  713. }
  714.         sfree(buffer);
  715. }
  717. seqloc_num = 0;
  718. if (program_number != blast_type_blastn)
  719. {
  720. if (do_default || do_seg)
  721. {
  722. SeqBufferSeg(sequence, length, offset, sparamsp, &seg_loc);
  723. SegParametersFree(sparamsp);
  724. sparamsp = NULL;
  725. seqloc_num++;
  726. }
  727. #ifdef CC_FILTER_ALLOWED
  728. if (do_coil_coil)
  729. {
  730. pccp = PccDatNew ();
  731. pccp->window = window_cc;
  732. ReadPccData (pccp);
  733. scores = PredictCCSeqLoc(slp, pccp);
  734. cc_slp = FilterCC(scores, cutoff_cc, length, linker_cc,
  735.                                           SeqIdDup(sip), FALSE);
  736. sfree(scores);
  737. PccDatFree (pccp);
  738. seqloc_num++;
  739. }
  740. #endif
  741. }
  742. else
  743. {
  744. if (do_default || do_dust)
  745. {
  746.          SeqBufferDust(sequence, length, offset, level_dust, window_dust,
  747.                        minwin_dust, linker_dust, &dust_loc);
  748. seqloc_num++;
  749. }
  750. #ifdef TEMP_BLAST_OPTIONS
  751. if (do_repeats)
  752. {
  753. /* Either the SeqLocPtr is SEQLOC_WHOLE (both strands) or SEQLOC_INT (probably 
  754. one strand).  In that case we make up a double-stranded one as we wish to look at both strands. */
  755. myslp_allocated = FALSE;
  756. if (slp->choice == SEQLOC_INT)
  757. {
  758. myslp = SeqLocIntNew(SeqLocStart(slp), SeqLocStop(slp), Seq_strand_both, SeqLocId(slp));
  759. myslp_allocated = TRUE;
  760. }
  761. else
  762. {
  763. myslp = slp;
  764. }
  765. start_timer;
  766. repeat_slp = BioseqHitRangeEngineByLoc(myslp, "blastn", repeat_database, repeat_options, NULL, NULL, NULL, NULL, NULL, 0);
  767. stop_timer("after repeat filtering");
  768. repeat_options = BLASTOptionDelete(repeat_options);
  769. sfree(repeat_database);
  770. if (myslp_allocated)
  771. SeqLocFree(myslp);
  772. seqloc_num++;
  773. }
  774. if (do_vecscreen)
  775. {
  776. /* Either the SeqLocPtr is SEQLOC_WHOLE (both strands) or SEQLOC_INT (probably 
  777. one strand).  In that case we make up a double-stranded one as we wish to look at both strands. */
  778. myslp_allocated = FALSE;
  779. if (slp->choice == SEQLOC_INT)
  780. {
  781. myslp = SeqLocIntNew(SeqLocStart(slp), SeqLocStop(slp), Seq_strand_both, SeqLocId(slp));
  782. myslp_allocated = TRUE;
  783. }
  784. else
  785. {
  786. myslp = slp;
  787. }
  788. VSScreenSequenceByLoc(myslp, NULL, vs_database, &seqalign, &vnp, NULL, NULL);
  789. vnp_var = vnp;
  790. while (vnp_var)
  791. {
  792. seqloc_tmp = vnp_var->ptr;
  793. if (vs_slp == NULL)
  794. {
  795. vs_slp = seqloc_tmp;
  796. }
  797. else
  798. {
  799. seqloc_var = vs_slp;
  800. while (seqloc_var->next)
  801. seqloc_var = seqloc_var->next;
  802. seqloc_var->next = seqloc_tmp;
  803. }
  804. vnp_var->ptr = NULL;
  805. vnp_var = vnp_var->next;
  806. }
  807. vnp = ListNodeFree(vnp);
  808. seqalign = SeqAlignSetFree(seqalign);
  809. sfree(vs_database);
  810. if (myslp_allocated)
  811. SeqLocFree(myslp);
  812. seqloc_num++;
  813. }
  814. #endif
  815. }
  817. if (seqloc_num)
  818. { 
  819. BlastSeqLoc* seqloc_list=NULL;  /* Holds all SeqLoc's for
  820.                                                       return. */
  821. #if 0
  822. if (seg_slp)
  823. ListNodeAddPointer(&seqloc_list, SEQLOC_MIX, seg_slp);
  824. if (cc_slp)
  825. ListNodeAddPointer(&seqloc_list, SEQLOC_MIX, cc_slp);
  826. if (dust_slp)
  827. ListNodeAddPointer(&seqloc_list, SEQLOC_MIX, dust_slp);
  828. if (repeat_slp)
  829. ListNodeAddPointer(&seqloc_list, SEQLOC_MIX, repeat_slp);
  830. if (vs_slp)
  831. ListNodeAddPointer(&seqloc_list, SEQLOC_MIX, vs_slp);
  832. #endif
  833.       if (dust_loc)
  834.          seqloc_list = dust_loc;
  835.       if (seg_loc)
  836.          seqloc_list = seg_loc;
  838. *seqloc_retval = seqloc_list;
  839. }
  841. return status;
  842. }
  844. static Int2
  845. GetFilteringLocationsForOneContext(BLAST_SequenceBlk* query_blk, BlastQueryInfo* query_info, Int2 context, Uint1 program_number, const char* filter_string, BlastSeqLoc* *filter_out, Boolean* mask_at_hash)
  846. {
  847.         Int2 status = 0;
  848.         Int4 query_length = 0;      /* Length of query described by SeqLocPtr. */
  849.         Int4 context_offset;
  850.         BlastMaskLoc *mask_slp, *mask_slp_var; /* Auxiliary locations for lower-case masking  */
  851.         BlastSeqLoc *filter_slp = NULL;     /* SeqLocPtr computed for filtering. */
  852.         BlastSeqLoc *filter_slp_combined;   /* Used to hold combined SeqLoc's */
  853.         Uint1 *buffer;              /* holds sequence for plus strand or protein. */
  855.         Boolean is_na = (program_number == blast_type_blastn);
  856.         Int2 index = (is_na ? context / 2 : context);
  858.         context_offset = query_info->context_offsets[context];
  859.         buffer = &query_blk->sequence[context_offset];
  861.         if ((query_length = BLAST_GetQueryLength(query_info, context)) <= 0)
  862.            return 0;
  865.         if ((status = BlastSetUp_Filter(program_number, buffer,
  866.                        query_length, 0, filter_string,
  867.                              mask_at_hash, &filter_slp))) 
  868.              return status;
  870.         /* Extract the mask locations corresponding to this query 
  871.                (frame, strand), detach it from other masks.
  872.                NB: for translated search the mask locations are expected in 
  873.                protein coordinates. The nucleotide locations must be converted
  874.                to protein coordinates prior to the call to BLAST_MainSetUp.
  875.         */
  876.         mask_slp = NULL;
  877.         for (mask_slp_var=query_blk->lcase_mask; mask_slp_var; mask_slp_var=mask_slp_var->next)
  878.         {
  879.                 if (mask_slp_var->index == index)
  880.                 {
  881.                    mask_slp = mask_slp_var;
  882.                    break;
  883.                 }
  884.         }
  886.         /* Attach the lower case mask locations to the filter locations and combine them */
  887.         if (mask_slp) {
  888.              if (filter_slp) {
  889.                   BlastSeqLoc *loc;           /* Iterator variable */
  890.                   for (loc = filter_slp; loc->next; loc = loc->next);
  891.                     loc->next = mask_slp->loc_list;
  892.              } else {
  893.                    filter_slp = mask_slp->loc_list;
  894.              }
  895.                 /* Set location list to NULL, to allow safe memory deallocation */
  896.                 mask_slp->loc_list = NULL;
  897.         }
  899.         filter_slp_combined = NULL;
  900.         CombineMaskLocations(filter_slp, &filter_slp_combined, 0);
  901.         *filter_out = filter_slp_combined;
  903.         filter_slp = BlastSeqLocFree(filter_slp);
  905. return 0;
  906. }
  909. Int2
  910. BlastSetUp_GetFilteringLocations(BLAST_SequenceBlk* query_blk, BlastQueryInfo* query_info, Uint1 program_number, const char* filter_string, BlastMaskLoc* *filter_out, Boolean* mask_at_hash, Blast_Message * *blast_message)
  911. {
  913.     Int2 status = 0;
  914.     Int2 context = 0; /* loop variable. */
  915.     const Boolean k_is_na = (program_number == blast_type_blastn);
  916.     BlastMaskLoc *last_maskloc = NULL;
  917.     BlastMaskLoc *filter_maskloc = NULL;   /* Local variable for mask locs. */
  918.     Boolean no_forward_strand = (query_info->first_context > 0);  /* filtering needed on reverse strand. */
  920.     for (context = query_info->first_context;
  921.          context <= query_info->last_context; ++context) {
  922.       
  923.         BlastSeqLoc *filter_per_context = NULL;   /* Used to hold combined SeqLoc's */
  924.         Boolean reverse = (k_is_na && ((context & 1) != 0));
  925.         Int4 query_length;
  927.         /* For each query, check if forward strand is present */
  928.         if ((query_length = BLAST_GetQueryLength(query_info, context)) < 0)
  929.         {
  930.             if ((context & 1) == 0)
  931.                no_forward_strand = TRUE;
  932.             continue;
  933.         }
  935.       if (!reverse || no_forward_strand)
  936.       {
  937.         if ((status=GetFilteringLocationsForOneContext(query_blk, query_info, context, program_number, filter_string, &filter_per_context, mask_at_hash)))
  938.         {
  939.                Blast_MessageWrite(blast_message, BLAST_SEV_ERROR, 2, 1, 
  940.                   "Failure at filtering");
  941.                return status;
  942.         }
  944.         /* NB: for translated searches filter locations are returned in 
  945.                protein coordinates, because the DNA lengths of sequences are 
  946.                not available here. The caller must take care of converting 
  947.                them back to nucleotide coordinates. */
  948.        if (filter_per_context)
  949.        {
  950.         if (!last_maskloc) {
  951.                 last_maskloc = (BlastMaskLoc *) calloc(1, sizeof(BlastMaskLoc));
  952.                 filter_maskloc = last_maskloc;
  953.         } else {
  954.                 last_maskloc->next =
  955.                         (BlastMaskLoc *) calloc(1, sizeof(BlastMaskLoc));
  956.                 last_maskloc = last_maskloc->next;
  957.         }
  959.         last_maskloc->index = (k_is_na ? context / 2 : context);
  960.         last_maskloc->loc_list = filter_per_context;
  961.        }
  962.       }
  963.     }
  965.     if (filter_out && filter_maskloc)
  966. *filter_out = filter_maskloc;
  968.     return 0;
  969. }
  971. Int2
  972. Blast_MaskTheResidues(Uint1 * buffer, Int4 length, Boolean is_na,
  973.                            ListNode * mask_loc, Boolean reverse, Int4 offset)
  974. {
  975.     SSeqRange *loc = NULL;
  976.     Int2 status = 0;
  977.     Int4 index, start, stop;
  978.     Uint1 mask_letter;
  980.     if (is_na)
  981.         mask_letter = 14;
  982.     else
  983.         mask_letter = 21;
  985.     for (; mask_loc; mask_loc = mask_loc->next) {
  986.         loc = (SSeqRange *) mask_loc->ptr;
  987.         if (reverse) {
  988.             start = length - 1 - loc->right;
  989.             stop = length - 1 - loc->left;
  990.         } else {
  991.             start = loc->left;
  992.             stop = loc->right;
  993.         }
  995.         start -= offset;
  996.         stop -= offset;
  998.         for (index = start; index <= stop; index++)
  999.             buffer[index] = mask_letter;
  1000.     }
  1002.     return status;
  1003. }
  1005. Int2 
  1006. BlastSetUp_MaskQuery(BLAST_SequenceBlk* query_blk, BlastQueryInfo* query_info, BlastMaskLoc *filter_maskloc, Uint1 program_number)
  1007. {
  1008.     const Boolean k_is_na = (program_number == blast_type_blastn);
  1009.     Int4 context; /* loop variable. */
  1010.     Int2 status=0;
  1012.     for (context = query_info->first_context;
  1013.          context <= query_info->last_context; ++context) {
  1014.       
  1015.         BlastMaskLoc* filter_maskloc_var = NULL;
  1016.         BlastSeqLoc *filter_per_context = NULL;   /* Used to hold combined SeqLoc's */
  1017.         Boolean reverse = (k_is_na && ((context & 1) != 0));
  1018.         Int4 query_length;
  1019.         Int4 context_offset;
  1020.         Uint1 *buffer;              /* holds sequence */
  1022.         /* For each query, check if forward strand is present */
  1023.         if ((query_length = BLAST_GetQueryLength(query_info, context)) < 0)
  1024.             continue;
  1026.         context_offset = query_info->context_offsets[context];
  1027.         buffer = &query_blk->sequence[context_offset];
  1029. filter_maskloc_var = filter_maskloc;
  1030.         while (filter_maskloc_var)
  1031.         {
  1032.              if (filter_maskloc_var->index == (k_is_na ? context / 2 : context))
  1033.              {
  1034. filter_per_context = filter_maskloc_var->loc_list;
  1035.                 break;
  1036.              }
  1037.              filter_maskloc_var = filter_maskloc_var->next;
  1038.         }
  1040.         if (buffer) {
  1042.             if ((status =
  1043.                      Blast_MaskTheResidues(buffer, query_length, k_is_na,
  1044.                                                 filter_per_context, reverse, 0)))
  1045.             {
  1046.                     return status;
  1047.             }
  1048.         }
  1049.     }
  1051.     return 0;
  1052. }