blast_inline.h
上传用户:yhdzpy8989
上传日期:2007-06-13
资源大小:13604k
文件大小:12k
源码类别:

生物技术

开发平台:

C/C++

  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: blast_inline.h,v $
  4.  * PRODUCTION Revision 1000.0  2004/06/01 18:13:43  gouriano
  5.  * PRODUCTION PRODUCTION: IMPORTED [GCC34_MSVC7] Dev-tree R1.2
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  9. /* $Id: blast_inline.h,v 1000.0 2004/06/01 18:13:43 gouriano Exp $
  10.  * ===========================================================================
  11.  *
  12.  *                            PUBLIC DOMAIN NOTICE
  13.  *               National Center for Biotechnology Information
  14.  *
  15.  *  This software/database is a "United States Government Work" under the
  16.  *  terms of the United States Copyright Act.  It was written as part of
  17.  *  the author's offical duties as a United States Government employee and
  18.  *  thus cannot be copyrighted.  This software/database is freely available
  19.  *  to the public for use. The National Library of Medicine and the U.S.
  20.  *  Government have not placed any restriction on its use or reproduction.
  21.  *
  22.  *  Although all reasonable efforts have been taken to ensure the accuracy
  23.  *  and reliability of the software and data, the NLM and the U.S.
  24.  *  Government do not and cannot warrant the performance or results that
  25.  *  may be obtained by using this software or data. The NLM and the U.S.
  26.  *  Government disclaim all warranties, express or implied, including
  27.  *  warranties of performance, merchantability or fitness for any particular
  28.  *  purpose.
  29.  *
  30.  *  Please cite the author in any work or product based on this material.
  31.  *
  32.  * ===========================================================================
  33.  *
  34.  */
  35. /** @file blast_inline.h
  36.  * @todo FIXME needs file description
  37.  */
  38. #include <algo/blast/core/mb_lookup.h>
  39. #include <algo/blast/core/blast_util.h>
  40. /** Given a word packed into an integer, compute a discontiguous word lookup 
  41.  *  index.
  42.  * @param subject Pointer to the next byte of the sequence after the end of 
  43.  *        the word (needed when word template is longer than 16 bases) [in]
  44.  * @param word A piece of the sequence packed into an integer [in]
  45.  * @param template_type What type of discontiguous word template to use [in]
  46.  * @return The lookup table index of the discontiguous word [out]
  47.  */
  48. static NCBI_INLINE Int4 ComputeDiscontiguousIndex(Uint1* subject, Int4 word,
  49.                   Uint1 template_type)
  50. {
  51.    Int4 index;
  52.    Int4 extra_code;   
  53.    switch (template_type) {
  54.    case TEMPL_11_16:
  55.       index = GET_WORD_INDEX_11_16(word);
  56.       break;
  57.    case TEMPL_12_16:
  58.       index = GET_WORD_INDEX_12_16(word);
  59.       break;
  60.    case TEMPL_11_16_OPT:
  61.       index = GET_WORD_INDEX_11_16_OPT(word);
  62.       break;
  63.    case TEMPL_12_16_OPT:
  64.       index = GET_WORD_INDEX_12_16_OPT(word);
  65.       break;
  66.    case TEMPL_11_18: 
  67.      extra_code = (Int4) GET_EXTRA_CODE_PACKED_4_18(subject);
  68.      index = (GET_WORD_INDEX_11_18(word) | extra_code);
  69.      break;
  70.    case TEMPL_12_18: 
  71.       extra_code = (Int4) GET_EXTRA_CODE_PACKED_4_18(subject);
  72.       index = (GET_WORD_INDEX_12_18(word) | extra_code);
  73.       break;
  74.    case TEMPL_11_18_OPT: 
  75.       extra_code = (Int4) GET_EXTRA_CODE_PACKED_4_18_OPT(subject);
  76.       index = (GET_WORD_INDEX_11_18_OPT(word) | extra_code);
  77.       break;
  78.    case TEMPL_12_18_OPT:
  79.       extra_code = (Int4) GET_EXTRA_CODE_PACKED_4_18_OPT(subject);
  80.       index = (GET_WORD_INDEX_12_18_OPT(word) | extra_code);
  81.       break;
  82.    case TEMPL_11_21: 
  83.       extra_code = (Int4) GET_EXTRA_CODE_PACKED_4_21(subject);
  84.       index = (GET_WORD_INDEX_11_21(word) | extra_code);
  85.       break;
  86.    case TEMPL_12_21:
  87.       extra_code = (Int4) GET_EXTRA_CODE_PACKED_4_21(subject);
  88.       index = (GET_WORD_INDEX_12_21(word) | extra_code);
  89.       break;
  90.    case TEMPL_11_21_OPT: 
  91.       extra_code = (Int4) GET_EXTRA_CODE_PACKED_4_21_OPT(subject);
  92.       index = (GET_WORD_INDEX_11_21_OPT(word) | extra_code);
  93.       break;
  94.    case TEMPL_12_21_OPT:
  95.       extra_code = (Int4) GET_EXTRA_CODE_PACKED_4_21_OPT(subject);
  96.       index = (GET_WORD_INDEX_12_21_OPT(word) | extra_code);
  97.          break;
  98.    default: 
  99.       extra_code = 0; 
  100.       index = 0;
  101.       break;
  102.    }
  103. #ifdef USE_HASH_TABLE
  104.    hash_buf = (Uint1*)&index;
  105.    CRC32(crc, hash_buf);
  106.    index = (crc>>hash_shift) & hash_mask;
  107. #endif
  108.    return index;
  109. }
  110. /** Compute the lookup table index for the first word template, given a word 
  111.  * position, template type and previous value of the word, in case of 
  112.  * one-base (2 bit) database scanning.
  113.  * @param word_start Pointer to the start of a word in the sequence [in]
  114.  * @param word The word packed into an integer value [in]
  115.  * @param sequence_bit By how many bits the real word start is shifted within 
  116.  *        a compressed sequence byte [in]
  117.  * @param template_type What discontiguous word template to use for index 
  118.  *        computation [in]
  119.  * @return The lookup index for the discontiguous word.
  120. */
  121. static NCBI_INLINE Int4 ComputeDiscontiguousIndex_1b(const Uint1* word_start, 
  122.                       Int4 word, Uint1 sequence_bit, Uint1 template_type)
  123. {
  124.    Int4 index;
  125.    Uint1* subject = (Uint1 *) word_start;
  126.    Uint1 bit;
  127.    Int4 extra_code, tmpval;   
  128.    /* Prepare auxiliary variables for extra code calculation */
  129.    tmpval = 0;
  130.    extra_code = 0;
  131.    /* The bits in an integer byte are counted in a reverse order than in a
  132.       sequence byte */
  133.    bit = 6 - sequence_bit;
  134.    switch (template_type) {
  135.    case TEMPL_11_16:
  136.       index = GET_WORD_INDEX_11_16(word);
  137.       break;
  138.    case TEMPL_12_16:
  139.       index = GET_WORD_INDEX_12_16(word);
  140.       break;
  141.    case TEMPL_11_16_OPT:
  142.       index = GET_WORD_INDEX_11_16_OPT(word);
  143.       break;
  144.    case TEMPL_12_16_OPT:
  145.       index = GET_WORD_INDEX_12_16_OPT(word);
  146.       break;
  147.    case TEMPL_11_18: 
  148.       GET_EXTRA_CODE_PACKED_18(subject, bit, tmpval, extra_code);
  149.       index = (GET_WORD_INDEX_11_18(word) | extra_code);
  150.       break;
  151.    case TEMPL_12_18: 
  152.       GET_EXTRA_CODE_PACKED_18(subject, bit, tmpval, extra_code);
  153.       index = (GET_WORD_INDEX_12_18(word) | extra_code);
  154.       break;
  155.    case TEMPL_11_18_OPT: 
  156.       GET_EXTRA_CODE_PACKED_18_OPT(subject, bit, tmpval, extra_code);
  157.       index = (GET_WORD_INDEX_11_18_OPT(word) | extra_code);
  158.       break;
  159.    case TEMPL_12_18_OPT:
  160.       GET_EXTRA_CODE_PACKED_18_OPT(subject, bit, tmpval, extra_code);
  161.       index = (GET_WORD_INDEX_12_18_OPT(word) | extra_code);
  162.       break;
  163.    case TEMPL_11_21: 
  164.       GET_EXTRA_CODE_PACKED_21(subject, bit, tmpval, extra_code);
  165.       index = (GET_WORD_INDEX_11_21(word) | extra_code);
  166.       break;
  167.    case TEMPL_12_21:
  168.       GET_EXTRA_CODE_PACKED_21(subject, bit, tmpval, extra_code);
  169.       index = (GET_WORD_INDEX_12_21(word) | extra_code);
  170.       break;
  171.    case TEMPL_11_21_OPT: 
  172.       GET_EXTRA_CODE_PACKED_21_OPT(subject, bit, tmpval, extra_code);
  173.       index = (GET_WORD_INDEX_11_21_OPT(word) | extra_code);
  174.       break;
  175.    case TEMPL_12_21_OPT:
  176.       GET_EXTRA_CODE_PACKED_21_OPT(subject, bit, tmpval, extra_code);
  177.       index = (GET_WORD_INDEX_12_21_OPT(word) | extra_code);
  178.          break;
  179.    default: 
  180.       extra_code = 0; 
  181.       index = 0;
  182.       break;
  183.    }
  184. #ifdef USE_HASH_TABLE
  185.    hash_buf = (Uint1*)&index;
  186.    CRC32(crc, hash_buf);
  187.    index = (crc>>hash_shift) & hash_mask;
  188. #endif
  189.   
  190.    return index;
  191. }
  192. static NCBI_INLINE void  _ComputeIndex(Int4 wordsize,
  193.   Int4 charsize,
  194.   Int4 mask,
  195.   const Uint1* word,
  196.   Int4* index);
  197. static NCBI_INLINE void  _ComputeIndexIncremental(Int4 wordsize,
  198.      Int4 charsize,
  199.      Int4 mask,
  200.      const Uint1* word,
  201.      Int4* index);
  202. /** Given a word, compute its index value from scratch.
  203.  *
  204.  * @param wordsize length of the word, in residues [in]
  205.  * @param charsize length of one residue, in bits [in]
  206.  * @param mask value used to mask the index so that only the bottom wordsize * charsize bits remain [in]
  207.  * @param word pointer to the beginning of the word [in]
  208.  * @param index the computed index value [out] 
  209.  */
  210. static NCBI_INLINE void  _ComputeIndex(Int4 wordsize,
  211.   Int4 charsize,
  212.   Int4 mask,
  213.   const Uint1* word,
  214.   Int4* index)
  215. {
  216.   Int4 i;
  217.   *index = 0;
  218.   for(i=0;i<wordsize;i++)
  219. {
  220.     *index = ((*index << charsize) | word[i]) & mask;
  221. }
  222.   return;
  223. }
  224. /** Given a word, compute its index value, reusing a previously 
  225.  *  computed index value.
  226.  *
  227.  * @param wordsize length of the word - 1, in residues [in]
  228.  * @param charsize length of one residue, in bits [in]
  229.  * @param mask value used to mask the index so that only the bottom wordsize * charsize bits remain [in]
  230.  * @param word pointer to the beginning of the word [in]
  231.  * @param index the computed index value [in/out]
  232.  */
  233. static NCBI_INLINE void  _ComputeIndexIncremental(Int4 wordsize,
  234.      Int4 charsize,
  235.      Int4 mask,
  236.      const Uint1* word,
  237.      Int4* index)
  238. {
  239.   *index = ((*index << charsize) | word[wordsize - 1]) & mask;
  240.   return;
  241. }
  242. /* Given a starting position of a word in a compressed nucleotide sequence, 
  243.  * compute this word's lookup table index
  244.  */
  245. static NCBI_INLINE Uint1* BlastNaLookupInitIndex(Int4 length,
  246.           const Uint1* word, Int4* index)
  247. {
  248.    Int4 i;
  249.    
  250.    *index = 0;
  251.    for (i = 0; i < length; ++i)
  252.       *index = ((*index)<<FULL_BYTE_SHIFT) | word[i];
  253.    return (Uint1 *) (word + length);
  254. }
  255. /* Recompute the word index given its previous value and the new location 
  256.  *  of the last byte of the word
  257.  */
  258. static NCBI_INLINE Int4 BlastNaLookupComputeIndex(Int4 scan_shift, Int4 mask, 
  259.       const Uint1* word, Int4 index)
  260. {
  261.    return (((index)<<scan_shift) & mask) | *(word);  
  262.    
  263. }
  264. /* Given a word computed from full bytes of a compressed sequence, 
  265.  *  shift it by 0-3 bases 
  266.  */
  267. static NCBI_INLINE Int4 BlastNaLookupAdjustIndex(Uint1* s, Int4 index, 
  268.                       Int4 mask, Uint1 bit)
  269. {
  270.    return (((index)<<bit) & mask) | ((*s)>>(FULL_BYTE_SHIFT-bit));
  271. }
  272. #define BLAST2NA_MASK 0xfc
  273. /** Compute the lookup index for a word in an uncompressed sequence, without 
  274.  *  using any previous index information.
  275.  * @param lookup Pointer to the traditional BLASTn lookup table structure [in]
  276.  * @param word Pointer to the start of the word [in]
  277.  * @param index The lookup index [out]
  278.  */
  279. static NCBI_INLINE Int2
  280. Na_LookupComputeIndex(LookupTable* lookup, Uint1* word, Int4* index)
  281. {
  282.    Int4 i;
  283.    Int4 wordsize = lookup->reduced_wordsize*COMPRESSION_RATIO; /* i.e. 8 or 4 */
  284.    *index = 0;
  285.    for (i = 0; i < wordsize; ++i) {
  286.       if ((word[i] & BLAST2NA_MASK) != 0) {
  287.  *index = 0;
  288.  return -1;
  289.       } else {
  290.  *index = (((*index)<<lookup->charsize) & lookup->mask) | word[i];
  291.       }
  292.    }
  293.    return 0;
  294. }
  295. /** Pack 4 sequence bytes into a one byte integer, assuming sequence contains
  296.  *  no ambiguities.
  297.  */
  298. #define PACK_WORD(q) ((q[0]<<6) + (q[1]<<4) + (q[2]<<2) + q[3]) 
  299. /** Compare a given number of bytes of an compressed subject sequence with
  300.  *  the non-compressed query sequence.
  301.  * @param q Pointer to the first byte to be compared in the query sequence [in]
  302.  * @param s Pointer to the first byte to be compared in the subject 
  303.  *        sequence [in]
  304.  * @param extra_bytes Number of compressed bytes to compare [in]
  305.  * @return TRUE if sequences are identical, FALSE if mismatch is found.
  306.  */
  307. static NCBI_INLINE Boolean BlastNaCompareExtraBytes(Uint1* q, Uint1* s, 
  308.                                                Int4 extra_bytes)
  309. {
  310.    Int4 index;
  311.    
  312.    for (index = 0; index < extra_bytes; ++index) {
  313.       if (*s++ != PACK_WORD(q))
  314.  return FALSE;
  315.       q += COMPRESSION_RATIO;
  316.    }
  317.    return TRUE;
  318. }
  319. /** Perform mini extension (up to max_left <= 4 bases) to the left; 
  320.  * @param q Pointer to the query base right after the ones to be extended [in]
  321.  * @param s Pointer to a byte in the compressed subject sequence that is to be 
  322.  *        tested for extension [in]
  323.  * @param max_left Maximal number of bits to compare [in]
  324.  * @return Number of matched bases
  325. */
  326. static NCBI_INLINE Uint1 
  327. BlastNaMiniExtendLeft(Uint1* q, const Uint1* s, Uint1 max_left)
  328. {
  329.    Uint1 left = 0;
  330.    for (left = 0; left < max_left; ++left) {
  331.       if (NCBI2NA_UNPACK_BASE(*s, left) != *--q) {
  332.  break;
  333.       }
  334.    }
  335.    return left;
  336. }
  337. /** Perform mini extension (up to max_right <= 4 bases) to the right; 
  338.  * @param q Pointer to the start of the extension in the query [in]
  339.  * @param s Pointer to a byte in the compressed subject sequence that is to be 
  340.  *        tested for extension [in]
  341.  * @param max_right Maximal number of bits to compare [in]
  342.  * @return Number of matched bases
  343. */
  344. static NCBI_INLINE Uint1 
  345. BlastNaMiniExtendRight(Uint1* q, const Uint1* s, Uint1 max_right)
  346. {
  347.    Uint1 right;
  348.    Uint1 index = 3;
  349.    
  350.    for (right = 0; right < max_right; ++right, --index) {
  351.       if (NCBI2NA_UNPACK_BASE(*s, index) != *q++) {
  352.  break;
  353.       }
  354.    }
  355.    return right;
  356. }