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

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

生物技术

开发平台：

C/C++

blast_setup_cxx.cpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: blast_setup_cxx.cpp,v $
  4.  * PRODUCTION Revision 1000.7  2004/06/01 18:05:58  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.68
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /*  $Id: blast_setup_cxx.cpp,v 1000.7 2004/06/01 18:05:58 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 official 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:  Christiam Camacho
  36. *
  37. * ===========================================================================
  38. */
  40. /// @file blast_setup_cxx.cpp
  41. /// Auxiliary setup functions for Blast objects interface.
  43. #include <ncbi_pch.hpp>
  44. #include <corelib/ncbiapp.hpp>
  45. #include <corelib/ncbireg.hpp>
  46. #include <corelib/metareg.hpp>
  47. #include <corelib/ncbifile.hpp>
  49. #include <objmgr/object_manager.hpp>
  50. #include <objmgr/scope.hpp>
  51. #include <objmgr/seq_vector.hpp>
  52. #include <objmgr/util/sequence.hpp>
  54. #include <objects/seqloc/Seq_loc.hpp>
  55. #include <objects/seq/seqport_util.hpp>
  56. #include <objects/seqfeat/Genetic_code_table.hpp>
  57. #include <objects/seq/Seq_data.hpp>
  58. #include <objects/seq/NCBIstdaa.hpp>
  60. #include <algo/blast/api/blast_options.hpp>
  61. #include "blast_setup.hpp"
  63. #include <algo/blast/core/blast_util.h>
  64. #include <algo/blast/core/blast_encoding.h>
  66. #include <algorithm>
  68. /** @addtogroup AlgoBlast
  69.  *
  70.  * @{
  71.  */
  73. BEGIN_NCBI_SCOPE
  74. BEGIN_SCOPE(blast)
  75. USING_SCOPE(ncbi::objects);
  77. /// Now allows query concatenation
  78. void
  79. SetupQueryInfo(const TSeqLocVector& queries, const CBlastOptions& options, 
  80.                BlastQueryInfo** qinfo)
  81. {
  82.     ASSERT(qinfo);
  84.     // Allocate and initialize the query info structure
  85.     if ( !((*qinfo) = (BlastQueryInfo*) calloc(1, sizeof(BlastQueryInfo)))) {
  86.         NCBI_THROW(CBlastException, eOutOfMemory, "Query info");
  87.     }
  89.     EProgram prog = options.GetProgram();
  90.     unsigned int nframes = GetNumberOfFrames(prog);
  91.     (*qinfo)->num_queries = static_cast<int>(queries.size());
  92.     (*qinfo)->first_context = 0;
  93.     (*qinfo)->last_context = (*qinfo)->num_queries * nframes - 1;
  95.     // Allocate the various arrays of the query info structure
  96.     int* context_offsets = NULL;
  97.     if ( !(context_offsets = (int*)
  98.            malloc(sizeof(int) * ((*qinfo)->last_context + 2)))) {
  99.         NCBI_THROW(CBlastException, eOutOfMemory, "Context offsets array");
  100.     }
  101.     if ( !((*qinfo)->eff_searchsp_array = 
  102.            (Int8*) calloc((*qinfo)->last_context + 1, sizeof(Int8)))) {
  103.         NCBI_THROW(CBlastException, eOutOfMemory, "Search space array");
  104.     }
  105.     if ( !((*qinfo)->length_adjustments = 
  106.            (int*) calloc((*qinfo)->last_context + 1, sizeof(int)))) {
  107.         NCBI_THROW(CBlastException, eOutOfMemory, "Length adjustments array");
  108.     }
  110.     (*qinfo)->context_offsets = context_offsets;
  111.     context_offsets[0] = 0;
  113.     bool is_na = (prog == eBlastn) ? true : false;
  114.     bool translate = 
  115.         ((prog == eBlastx) || (prog == eTblastx)) ? true : false;
  117.     // Adjust first context depending on the first query strand
  118.     // Unless the strand option is set to single strand, the actual
  119.     // CSeq_locs dictate which strand to examine during the search.
  120.     ENa_strand strand_opt = options.GetStrandOption();
  121.     ENa_strand strand;
  123.     if (is_na || translate) {
  124.         if (strand_opt == eNa_strand_both || 
  125.             strand_opt == eNa_strand_unknown) {
  126.             strand = sequence::GetStrand(*queries.front().seqloc, 
  127.                                          queries.front().scope);
  128.         } else {
  129.             strand = strand_opt;
  130.         }
  132.         if (strand == eNa_strand_minus) {
  133.             if (translate) {
  134.                 (*qinfo)->first_context = 3;
  135.             } else {
  136.                 (*qinfo)->first_context = 1;
  137.             }
  138.         }
  139.     }
  141.     // Set up the context offsets into the sequence that will be added
  142.     // to the sequence block structure.
  143.     unsigned int ctx_index = 0;      // index into context_offsets array
  144.     // Longest query length, to be saved in the query info structure
  145.     Uint4 max_length = 0;
  147.     ITERATE(TSeqLocVector, itr, queries) {
  148.         TSeqPos length = sequence::GetLength(*itr->seqloc, itr->scope);
  149.         ASSERT(length != numeric_limits<TSeqPos>::max());
  151.         strand = sequence::GetStrand(*itr->seqloc, itr->scope);
  152.         if (strand_opt == eNa_strand_minus || strand_opt == eNa_strand_plus) {
  153.             strand = strand_opt;
  154.         }
  156.         if (translate) {
  157.             for (unsigned int i = 0; i < nframes; i++) {
  158.                 unsigned int prot_length = 
  159.                     (length - i % CODON_LENGTH) / CODON_LENGTH;
  160.                 max_length = MAX(max_length, prot_length);
  162.                 switch (strand) {
  163.                 case eNa_strand_plus:
  164.                     if (i < 3) {
  165.                         context_offsets[ctx_index + i + 1] = 
  166.                             context_offsets[ctx_index + i] + prot_length + 1;
  167.                     } else {
  168.                         context_offsets[ctx_index + i + 1] = 
  169.                             context_offsets[ctx_index + i];
  170.                     }
  171.                     break;
  173.                 case eNa_strand_minus:
  174.                     if (i < 3) {
  175.                         context_offsets[ctx_index + i + 1] = 
  176.                             context_offsets[ctx_index + i];
  177.                     } else {
  178.                         context_offsets[ctx_index + i + 1] =
  179.                             context_offsets[ctx_index + i] + prot_length + 1;
  180.                     }
  181.                     break;
  183.                 case eNa_strand_both:
  184.                 case eNa_strand_unknown:
  185.                     context_offsets[ctx_index + i + 1] = 
  186.                         context_offsets[ctx_index + i] + prot_length + 1;
  187.                     break;
  189.                 default:
  190.                     abort();
  191.                 }
  192.             }
  193.         } else {
  194.             max_length = MAX(max_length, length);
  195.             if (is_na) {
  196.                 switch (strand) {
  197.                 case eNa_strand_plus:
  198.                     context_offsets[ctx_index + 1] =
  199.                         context_offsets[ctx_index] + length + 1;
  200.                     context_offsets[ctx_index + 2] =
  201.                         context_offsets[ctx_index + 1];
  202.                     break;
  204.                 case eNa_strand_minus:
  205.                     context_offsets[ctx_index + 1] =
  206.                         context_offsets[ctx_index];
  207.                     context_offsets[ctx_index + 2] =
  208.                         context_offsets[ctx_index + 1] + length + 1;
  209.                     break;
  211.                 case eNa_strand_both:
  212.                 case eNa_strand_unknown:
  213.                     context_offsets[ctx_index + 1] =
  214.                         context_offsets[ctx_index] + length + 1;
  215.                     context_offsets[ctx_index + 2] =
  216.                         context_offsets[ctx_index + 1] + length + 1;
  217.                     break;
  219.                 default:
  220.                     abort();
  221.                 }
  222.             } else {    // protein
  223.                 context_offsets[ctx_index + 1] = 
  224.                     context_offsets[ctx_index] + length + 1;
  225.             }
  226.         }
  227.         ctx_index += nframes;
  228.     }
  229.     (*qinfo)->max_length = max_length;
  230. }
  232. Uint1
  233. GetQueryEncoding(EProgram program)
  234. {
  235.     Uint1 retval = 0;
  237.     switch (program) {
  238.     case eBlastn: 
  239.         retval = BLASTNA_ENCODING; 
  240.         break;
  242.     case eBlastp: 
  243.     case eTblastn:
  244.     case eRPSBlast: 
  245.         retval = BLASTP_ENCODING; 
  246.         break;
  248.     case eBlastx:
  249.     case eTblastx:
  250.     case eRPSTblastn:
  251.         retval = NCBI4NA_ENCODING;
  252.         break;
  254.     default:
  255.         NCBI_THROW(CBlastException, eBadParameter, "Query Encoding");
  256.     }
  258.     return retval;
  259. }
  261. Uint1
  262. GetSubjectEncoding(EProgram program)
  263. {
  264.     Uint1 retval = 0;
  266.     switch (program) {
  267.     case eBlastn: 
  268.         retval = BLASTNA_ENCODING; 
  269.         break;
  271.     case eBlastp: 
  272.     case eBlastx:
  273.         retval = BLASTP_ENCODING; 
  274.         break;
  276.     case eTblastn:
  277.     case eTblastx:
  278.         retval = NCBI4NA_ENCODING;
  279.         break;
  281.     default:
  282.         NCBI_THROW(CBlastException, eBadParameter, "Subject Encoding");
  283.     }
  285.     return retval;
  286. }
  288. void
  289. SetupQueries(const TSeqLocVector& queries, const CBlastOptions& options,
  290.              const CBlastQueryInfo& qinfo, BLAST_SequenceBlk** seqblk)
  291. {
  292.     ASSERT(seqblk);
  293.     ASSERT(queries.size() != 0);
  295.     EProgram prog = options.GetProgram();
  297.     // Determine sequence encoding
  298.     Uint1 encoding = GetQueryEncoding(prog);
  300.     int buflen = qinfo->context_offsets[qinfo->last_context+1] + 1;
  301.     Uint1* buf = (Uint1*) calloc((buflen+1), sizeof(Uint1));
  302.     if ( !buf ) {
  303.         NCBI_THROW(CBlastException, eOutOfMemory, "Query sequence buffer");
  304.     }
  306.     bool is_na = (prog == eBlastn) ? true : false;
  307.     bool translate = ((prog == eBlastx) || (prog == eTblastx)) ? true : false;
  309.     unsigned int ctx_index = 0;      // index into context_offsets array
  310.     unsigned int nframes = GetNumberOfFrames(prog);
  312.     BlastMaskLoc* mask = NULL, *head_mask = NULL, *last_mask = NULL;
  314.     // Unless the strand option is set to single strand, the actual
  315.     // CSeq_locs dictacte which strand to examine during the search
  316.     ENa_strand strand_opt = options.GetStrandOption();
  317.     int index = 0;
  319.     ITERATE(TSeqLocVector, itr, queries) {
  321.         ENa_strand strand;
  323.         if ((is_na || translate) &&
  324.             (strand_opt == eNa_strand_unknown || 
  325.              strand_opt == eNa_strand_both)) 
  326.         {
  327.             strand = sequence::GetStrand(*itr->seqloc, itr->scope);
  328.         } else {
  329.             strand = strand_opt;
  330.         }
  332.         if (itr->mask)
  333.             mask = CSeqLoc2BlastMaskLoc(*itr->mask, index);
  335.         ++index;
  337.         if (translate) {
  338.             ASSERT(strand == eNa_strand_both ||
  339.                    strand == eNa_strand_plus ||
  340.                    strand == eNa_strand_minus);
  342.             // Get both strands of the original nucleotide sequence with
  343.             // sentinels
  344.             pair<AutoPtr<Uint1, CDeleter<Uint1> >, TSeqPos> seqbuf(
  345.                 GetSequence(*itr->seqloc, encoding, itr->scope, strand,
  346.                             eSentinels));
  347.             
  349.             AutoPtr<Uint1, ArrayDeleter<Uint1> > gc = 
  350.                 FindGeneticCode(options.GetQueryGeneticCode());
  351.             int na_length = sequence::GetLength(*itr->seqloc, itr->scope);
  352.             Uint1* seqbuf_rev = NULL;  // negative strand
  353.             if (strand == eNa_strand_both)
  354.                seqbuf_rev = seqbuf.first.get() + na_length + 1;
  355.             else if (strand == eNa_strand_minus)
  356.                seqbuf_rev = seqbuf.first.get() + 1;
  358.             // Populate the sequence buffer
  359.             for (unsigned int i = 0; i < nframes; i++) {
  360.                 if (BLAST_GetQueryLength(qinfo, i) <= 0) {
  361.                     continue;
  362.                 }
  364.                 int offset = qinfo->context_offsets[ctx_index + i];
  365.                 short frame = BLAST_ContextToFrame(prog, i);
  366.                 BLAST_GetTranslation(seqbuf.first.get() + 1, seqbuf_rev,
  367.                    na_length, frame, &buf[offset], gc.get());
  368.             }
  369.             // Translate the lower case mask coordinates;
  370.             BlastMaskLocDNAToProtein(&mask, *itr->seqloc, itr->scope);
  372.         } else if (is_na) {
  374.             ASSERT(strand == eNa_strand_both ||
  375.                    strand == eNa_strand_plus ||
  376.                    strand == eNa_strand_minus);
  378.             pair<AutoPtr<Uint1, CDeleter<Uint1> >, TSeqPos> seqbuf(
  379.                 GetSequence(*itr->seqloc, encoding, itr->scope, strand,
  380.                             eSentinels));
  381.             int idx = (strand == eNa_strand_minus) ? 
  382.                 ctx_index + 1 : ctx_index;
  383.             int offset = qinfo->context_offsets[idx];
  384.             memcpy(&buf[offset], seqbuf.first.get(), seqbuf.second);
  386.         } else {
  388.             pair<AutoPtr<Uint1, CDeleter<Uint1> >, TSeqPos> seqbuf(
  389.                 GetSequence(*itr->seqloc, encoding, itr->scope,
  390.                             eNa_strand_unknown, eSentinels));
  391.             int offset = qinfo->context_offsets[ctx_index];
  392.             memcpy(&buf[offset], seqbuf.first.get(), seqbuf.second);
  394.         }
  396.         ctx_index += nframes;
  397.         
  398.         if (mask) {
  399.             if ( !last_mask )
  400.                 head_mask = last_mask = mask;
  401.             else {
  402.                 last_mask->next = mask;
  403.                 last_mask = last_mask->next;
  404.             }
  405.         }
  406.     }
  408.     if (BlastSeqBlkNew(seqblk) < 0) {
  409.         NCBI_THROW(CBlastException, eOutOfMemory, "Query sequence block");
  410.     }
  412.     BlastSeqBlkSetSequence(*seqblk, buf, buflen - 2);
  414.     (*seqblk)->lcase_mask = head_mask;
  415.     (*seqblk)->lcase_mask_allocated = TRUE;
  417.     return;
  418. }
  421. void
  422. SetupSubjects(const TSeqLocVector& subjects, 
  423.               EProgram prog,
  424.               vector<BLAST_SequenceBlk*>* seqblk_vec, 
  425.               unsigned int* max_subjlen)
  426. {
  427.     ASSERT(seqblk_vec);
  428.     ASSERT(max_subjlen);
  429.     ASSERT(subjects.size() != 0);
  431.     // Nucleotide subject sequences are stored in ncbi2na format, but the
  432.     // uncompressed format (ncbi4na/blastna) is also kept to re-evaluate with
  433.     // the ambiguities
  434.     bool subj_is_na = (prog == eBlastn  ||
  435.                        prog == eTblastn ||
  436.                        prog == eTblastx);
  438.     ESentinelType sentinels = eSentinels;
  439.     if (prog == eTblastn || prog == eTblastx) {
  440.         sentinels = eNoSentinels;
  441.     }
  443.     Uint1 encoding = GetSubjectEncoding(prog);
  444.        
  445.     // TODO: Should strand selection on the subject sequences be allowed?
  446.     //ENa_strand strand = options->GetStrandOption(); 
  447.     int index = 0; // Needed for lower case masks only.
  449.     *max_subjlen = 0;
  451.     ITERATE(TSeqLocVector, itr, subjects) {
  452.         BLAST_SequenceBlk* subj = NULL;
  454.         pair<AutoPtr<Uint1, CDeleter<Uint1> >, TSeqPos> seqbuf(
  455.             GetSequence(*itr->seqloc, encoding, itr->scope,
  456.                         eNa_strand_plus, sentinels));
  458.         if (BlastSeqBlkNew(&subj) < 0) {
  459.             NCBI_THROW(CBlastException, eOutOfMemory, "Subject sequence block");
  460.         }
  462.         /* Set the lower case mask, if it exists */
  463.         if (itr->mask)
  464.             subj->lcase_mask = CSeqLoc2BlastMaskLoc(*itr->mask, index);
  465.         ++index;
  467.         if (subj_is_na) {
  468.             BlastSeqBlkSetSequence(subj, seqbuf.first.release(), 
  469.                (sentinels == eSentinels) ? seqbuf.second - 2 : seqbuf.second);
  471.             try {
  472.                 // Get the compressed sequence
  473.                 pair<AutoPtr<Uint1, CDeleter<Uint1> >, TSeqPos> comp_seqbuf(
  474.                     GetSequence(*itr->seqloc, NCBI2NA_ENCODING, itr->scope,
  475.                                  eNa_strand_plus, eNoSentinels));
  476.                 BlastSeqBlkSetCompressedSequence(subj, 
  477.                                                  comp_seqbuf.first.release());
  478.             } catch (const CSeqVectorException&) {
  479.                 BlastSequenceBlkFree(subj);
  480.                 NCBI_THROW(CBlastException, eInvalidCharacter, 
  481.                            "Gaps found in subject sequence");
  482.             }
  483.         } else {
  484.             BlastSeqBlkSetSequence(subj, seqbuf.first.release(), 
  485.                                    seqbuf.second - 2);
  486.         }
  488.         seqblk_vec->push_back(subj);
  489.         (*max_subjlen) = MAX((*max_subjlen),
  490.                 sequence::GetLength(*itr->seqloc, itr->scope));
  492.     }
  493. }
  495. TSeqPos CalculateSeqBufferLength(TSeqPos sequence_length, Uint1 encoding,
  496.                                  ENa_strand strand, ESentinelType sentinel)
  497.                                  THROWS((CBlastException))
  498. {
  499.     TSeqPos retval = 0;
  501.     switch (encoding) {
  502.     // Strand and sentinels are irrelevant in this encoding.
  503.     // Strand is always plus and sentinels cannot be represented
  504.     case NCBI2NA_ENCODING:
  505.         ASSERT(sentinel == eNoSentinels);
  506.         ASSERT(strand == eNa_strand_plus);
  507.         retval = sequence_length / COMPRESSION_RATIO + 1;
  508.         break;
  510.     case NCBI4NA_ENCODING:
  511.         if (sentinel == eSentinels) {
  512.             if (strand == eNa_strand_both) {
  513.                 retval = sequence_length * 2;
  514.                 retval += 3;
  515.             } else {
  516.                 retval = sequence_length + 2;
  517.             }
  518.         } else {
  519.             if (strand == eNa_strand_both) {
  520.                 retval = sequence_length * 2;
  521.             } else {
  522.                 retval = sequence_length;
  523.             }
  524.         }
  525.         break;
  527.     case BLASTP_ENCODING:
  528.         ASSERT(sentinel == eSentinels);
  529.         ASSERT(strand == eNa_strand_unknown);
  530.         retval = sequence_length + 2;
  531.         break;
  533.     default:
  534.         NCBI_THROW(CBlastException, eBadParameter, "Unsupported encoding");
  535.     }
  537.     return retval;
  538. }
  540. // Compresses sequence data on vector to buffer, which should have been
  541. // allocated and have the right size.
  542. static 
  543. void CompressDNA(const CSeqVector& vec, Uint1* buffer, const int buflen)
  544. {
  545.     TSeqPos i;                  // loop index of original sequence
  546.     TSeqPos ci;                 // loop index for compressed sequence
  548.     ASSERT(vec.GetCoding() == CSeq_data::e_Ncbi2na);
  550.     for (ci = 0, i = 0; ci < (TSeqPos) buflen-1; ci++, i += COMPRESSION_RATIO) {
  551.         buffer[ci] = ((vec[i+0] & NCBI2NA_MASK)<<6) |
  552.                      ((vec[i+1] & NCBI2NA_MASK)<<4) |
  553.                      ((vec[i+2] & NCBI2NA_MASK)<<2) |
  554.                      ((vec[i+3] & NCBI2NA_MASK)<<0);
  555.     }
  557.     buffer[ci] = 0;
  558.     for (; i < vec.size(); i++) {
  559.             Uint1 bit_shift = 0;
  560.             switch (i%COMPRESSION_RATIO) {
  561.                case 0: bit_shift = 6; break;
  562.                case 1: bit_shift = 4; break;
  563.                case 2: bit_shift = 2; break;
  564.                default: abort();   // should never happen
  565.             }
  566.             buffer[ci] |= ((vec[i] & NCBI2NA_MASK)<<bit_shift);
  567.     }
  568.     // Set the number of bases in the last byte.
  569.     buffer[ci] |= vec.size()%COMPRESSION_RATIO;
  570. }
  572. Uint1 GetSentinelByte(Uint1 encoding) THROWS((CBlastException))
  573. {
  574.     switch (encoding) {
  575.     case BLASTP_ENCODING:
  576.         return NULLB;
  578.     case NCBI4NA_ENCODING:
  579.     case BLASTNA_ENCODING:
  580.         return 0xF;
  582.     default:
  583.         NCBI_THROW(CBlastException, eBadParameter, "Unsupported encoding");
  584.     }
  585. }
  587. pair<AutoPtr<Uint1, CDeleter<Uint1> >, TSeqPos>
  588. GetSequence(const CSeq_loc& sl, Uint1 encoding, CScope* scope,
  589.             ENa_strand strand, ESentinelType sentinel) 
  590.             THROWS((CBlastException, CException))
  591. {
  592.     Uint1* buf = NULL;          // buffer to write sequence
  593.     Uint1* buf_var = NULL;      // temporary pointer to buffer
  594.     TSeqPos buflen;             // length of buffer allocated
  595.     TSeqPos i;                  // loop index of original sequence
  596.     AutoPtr<Uint1, CDeleter<Uint1> > safe_buf; // contains buf to ensure 
  597.                                                // exception safety
  599.     CBioseq_Handle handle = scope->GetBioseqHandle(sl); // might throw exception
  601.     // Retrieves the correct strand (plus or minus), but not both
  602.     CSeqVector sv =
  603.         handle.GetSequenceView(sl, CBioseq_Handle::eViewConstructed,
  604.                                CBioseq_Handle::eCoding_Ncbi);
  606.     switch (encoding) {
  607.     // Protein sequences (query & subject) always have sentinels around sequence
  608.     case BLASTP_ENCODING:
  609.         sv.SetCoding(CSeq_data::e_Ncbistdaa);
  610.         buflen = CalculateSeqBufferLength(sv.size(), BLASTP_ENCODING);
  611.         buf = buf_var = (Uint1*) malloc(sizeof(Uint1)*buflen);
  612.         safe_buf.reset(buf);
  613.         *buf_var++ = GetSentinelByte(encoding);
  614.         for (i = 0; i < sv.size(); i++)
  615.             *buf_var++ = sv[i];
  616.         *buf_var++ = GetSentinelByte(encoding);
  617.         break;
  619.     case NCBI4NA_ENCODING:
  620.     case BLASTNA_ENCODING: // Used for nucleotide blastn queries
  621.         sv.SetCoding(CSeq_data::e_Ncbi4na);
  622.         buflen = CalculateSeqBufferLength(sv.size(), NCBI4NA_ENCODING,
  623.                                           strand, sentinel);
  624.         buf = buf_var = (Uint1*) malloc(sizeof(Uint1)*buflen);
  625.         safe_buf.reset(buf);
  626.         if (sentinel == eSentinels)
  627.             *buf_var++ = GetSentinelByte(encoding);
  629.         if (encoding == BLASTNA_ENCODING) {
  630.             for (i = 0; i < sv.size(); i++) {
  631.                 *buf_var++ = NCBI4NA_TO_BLASTNA[sv[i]];
  632.             }
  633.         } else {
  634.             for (i = 0; i < sv.size(); i++) {
  635.                 *buf_var++ = sv[i];
  636.             }
  637.         }
  638.         if (sentinel == eSentinels)
  639.             *buf_var++ = GetSentinelByte(encoding);
  641.         if (strand == eNa_strand_both) {
  642.             // Get the minus strand if both strands are required
  643.             sv = handle.GetSequenceView(sl, CBioseq_Handle::eViewConstructed,
  644.                                         CBioseq_Handle::eCoding_Ncbi, 
  645.                                         eNa_strand_minus);
  646.             sv.SetCoding(CSeq_data::e_Ncbi4na);
  647.             if (encoding == BLASTNA_ENCODING) {
  648.                 for (i = 0; i < sv.size(); i++) {
  649.                     *buf_var++ = NCBI4NA_TO_BLASTNA[sv[i]];
  650.                 }
  651.             } else {
  652.                 for (i = 0; i < sv.size(); i++) {
  653.                     *buf_var++ = sv[i];
  654.                 }
  655.             }
  656.             if (sentinel == eSentinels) {
  657.                 *buf_var++ = GetSentinelByte(encoding);
  658.             }
  659.         }
  661.         break;
  663.     /* Used only in Blast2Sequences for the subject sequence. 
  664.      * No sentinels can be used. As in readdb, remainder 
  665.      * (sv.size()%COMPRESSION_RATIO != 0) goes in the last 2 bits of the 
  666.      * last byte.
  667.      */
  668.     case NCBI2NA_ENCODING:
  669.         ASSERT(sentinel == eNoSentinels);
  670.         sv.SetCoding(CSeq_data::e_Ncbi2na);
  671.         buflen = CalculateSeqBufferLength(sv.size(), sv.GetCoding(),
  672.                                           eNa_strand_plus, eNoSentinels);
  673.         buf = (Uint1*) malloc(sizeof(Uint1)*buflen);
  674.         safe_buf.reset(buf);
  675.         CompressDNA(sv, buf, buflen);
  676.         break;
  678.     default:
  679.         NCBI_THROW(CBlastException, eBadParameter, "Invalid encoding");
  680.     }
  682.     return make_pair(safe_buf, buflen);
  683. }
  685. #if 0
  686. // Not used right now, need to complete implementation
  687. void
  688. BLASTGetTranslation(const Uint1* seq, const Uint1* seq_rev,
  689.         const int nucl_length, const short frame, Uint1* translation)
  690. {
  691.     TSeqPos ni = 0;     // index into nucleotide sequence
  692.     TSeqPos pi = 0;     // index into protein sequence
  694.     const Uint1* nucl_seq = frame >= 0 ? seq : seq_rev;
  695.     translation[0] = NULLB;
  696.     for (ni = ABS(frame)-1; ni < (TSeqPos) nucl_length-2; ni += CODON_LENGTH) {
  697.         Uint1 residue = CGen_code_table::CodonToIndex(nucl_seq[ni+0], 
  698.                                                       nucl_seq[ni+1],
  699.                                                       nucl_seq[ni+2]);
  700.         if (IS_residue(residue))
  701.             translation[pi++] = residue;
  702.     }
  703.     translation[pi++] = NULLB;
  705.     return;
  706. }
  707. #endif
  709. AutoPtr<Uint1, ArrayDeleter<Uint1> >
  710. FindGeneticCode(int genetic_code)
  711. {
  712.     Uint1* retval = NULL;
  713.     CSeq_data gc_ncbieaa(CGen_code_table::GetNcbieaa(genetic_code),
  714.             CSeq_data::e_Ncbieaa);
  715.     CSeq_data gc_ncbistdaa;
  717.     TSeqPos nconv = CSeqportUtil::Convert(gc_ncbieaa, &gc_ncbistdaa,
  718.             CSeq_data::e_Ncbistdaa);
  720.     ASSERT(gc_ncbistdaa.IsNcbistdaa());
  721.     ASSERT(nconv == gc_ncbistdaa.GetNcbistdaa().Get().size());
  723.     try {
  724.         retval = new Uint1[nconv];
  725.     } catch (const bad_alloc&) {
  726.         return NULL;
  727.     }
  729.     for (unsigned int i = 0; i < nconv; i++)
  730.         retval[i] = gc_ncbistdaa.GetNcbistdaa().Get()[i];
  732.     return retval;
  733. }
  735. string
  736. FindMatrixPath(const char* matrix_name, bool is_prot)
  737. {
  738.     string retval;
  739.     string full_path;       // full path to matrix file
  741.     if (!matrix_name)
  742.         return retval;
  744.     string mtx(matrix_name);
  745.     mtx = NStr::ToUpper(mtx);
  747.     // Look for matrix file in local directory
  748.     full_path = mtx;
  749.     if (CFile(full_path).Exists()) {
  750.         return retval;
  751.     }
  753.     // Obtain the matrix path from the ncbi configuration file
  754.     CMetaRegistry::SEntry sentry;
  755.     sentry = CMetaRegistry::Load("ncbi", CMetaRegistry::eName_RcOrIni);
  756.     string path = sentry.registry ? sentry.registry->Get("NCBI", "Data") : "";
  758.     full_path = CFile::MakePath(path, mtx);
  759.     if (CFile(full_path).Exists()) {
  760.         retval = full_path;
  761.         retval.erase(retval.size() - mtx.size());
  762.         return retval;
  763.     }
  765.     // Try appending "aa" or "nt" 
  766.     full_path = path;
  767.     full_path += CFile::AddTrailingPathSeparator(full_path);
  768.     full_path += is_prot ? "aa" : "nt";
  769.     full_path += CFile::AddTrailingPathSeparator(full_path);
  770.     full_path += mtx;
  771.     if (CFile(full_path).Exists()) {
  772.         retval = full_path;
  773.         retval.erase(retval.size() - mtx.size());
  774.         return retval;
  775.     }
  777.     // Try using local "data" directory
  778.     full_path = "data";
  779.     full_path += CFile::AddTrailingPathSeparator(full_path);
  780.     full_path += mtx;
  781.     if (CFile(full_path).Exists()) {
  782.         retval = full_path;
  783.         retval.erase(retval.size() - mtx.size());
  784.         return retval;
  785.     }
  787.     CNcbiApplication* app = CNcbiApplication::Instance();
  788.     if (!app)
  789.         return retval;
  791.     const string& blastmat_env = app->GetEnvironment().Get("BLASTMAT");
  792.     if (CFile(blastmat_env).Exists()) {
  793.         full_path = blastmat_env;
  794.         full_path += CFile::AddTrailingPathSeparator(full_path);
  795.         full_path += is_prot ? "aa" : "nt";
  796.         full_path += CFile::AddTrailingPathSeparator(full_path);
  797.         full_path += mtx;
  798.         if (CFile(full_path).Exists()) {
  799.             retval = full_path;
  800.             retval.erase(retval.size() - mtx.size());
  801.             return retval;
  802.         }
  803.     }
  805. #ifdef OS_UNIX
  806.     full_path = BLASTMAT_DIR;
  807.     full_path += CFile::AddTrailingPathSeparator(full_path);
  808.     full_path += is_prot ? "aa" : "nt";
  809.     full_path += CFile::AddTrailingPathSeparator(full_path);
  810.     full_path += mtx;
  811.     if (CFile(full_path).Exists()) {
  812.         retval = full_path;
  813.         retval.erase(retval.size() - mtx.size());
  814.         return retval;
  815.     }
  816. #endif
  818.     // Try again without the "aa" or "nt"
  819.     if (CFile(blastmat_env).Exists()) {
  820.         full_path = blastmat_env;
  821.         full_path += CFile::AddTrailingPathSeparator(full_path);
  822.         full_path += mtx;
  823.         if (CFile(full_path).Exists()) {
  824.             retval = full_path;
  825.             retval.erase(retval.size() - mtx.size());
  826.             return retval;
  827.         }
  828.     }
  830. #ifdef OS_UNIX
  831.     full_path = BLASTMAT_DIR;
  832.     full_path += CFile::AddTrailingPathSeparator(full_path);
  833.     full_path += mtx;
  834.     if (CFile(full_path).Exists()) {
  835.         retval = full_path;
  836.         retval.erase(retval.size() - mtx.size());
  837.         return retval;
  838.     }
  839. #endif
  841.     return retval;
  842. }
  844. /// Checks if a BLAST database exists at a given file path: looks for 
  845. /// an alias file first, then for an index file
  846. static bool BlastDbFileExists(string& path, bool is_prot)
  847. {
  848.     // Check for alias file first
  849.     string full_path = path + (is_prot ? ".pal" : ".nal");
  850.     if (CFile(full_path).Exists())
  851.         return true;
  852.     // Check for an index file
  853.     full_path = path + (is_prot ? ".pin" : ".nin");
  854.     if (CFile(full_path).Exists())
  855.         return true;
  856.     return false;
  857. }
  859. string
  860. FindBlastDbPath(const char* dbname, bool is_prot)
  861. {
  862.     string retval;
  863.     string full_path;       // full path to matrix file
  865.     if (!dbname)
  866.         return retval;
  868.     string database(dbname);
  870.     // Look for matrix file in local directory
  871.     full_path = database;
  872.     if (BlastDbFileExists(full_path, is_prot)) {
  873.         return retval;
  874.     }
  876.     CNcbiApplication* app = CNcbiApplication::Instance();
  877.     if (app) {
  878.         const string& blastdb_env = app->GetEnvironment().Get("BLASTDB");
  879.         if (CFile(blastdb_env).Exists()) {
  880.             full_path = blastdb_env;
  881.             full_path += CFile::AddTrailingPathSeparator(full_path);
  882.             full_path += database;
  883.             if (BlastDbFileExists(full_path, is_prot)) {
  884.                 retval = full_path;
  885.                 retval.erase(retval.size() - database.size());
  886.                 return retval;
  887.             }
  888.         }
  889.     }
  891.     // Obtain the matrix path from the ncbi configuration file
  892.     CMetaRegistry::SEntry sentry;
  893.     sentry = CMetaRegistry::Load("ncbi", CMetaRegistry::eName_RcOrIni);
  894.     string path = 
  895.         sentry.registry ? sentry.registry->Get("BLAST", "BLASTDB") : "";
  897.     full_path = CFile::MakePath(path, database);
  898.     if (BlastDbFileExists(full_path, is_prot)) {
  899.         retval = full_path;
  900.         retval.erase(retval.size() - database.size());
  901.         return retval;
  902.     }
  904.     return retval;
  905. }
  907. unsigned int
  908. GetNumberOfFrames(EProgram p)
  909. {
  910.     unsigned int retval = 0;
  912.     switch (p) {
  913.     case eBlastn:
  914.         retval = 2;
  915.         break;
  916.     case eBlastp:
  917.     case eRPSBlast:
  918.     case eTblastn: 
  919.     case eRPSTblastn: 
  920.         retval = 1;
  921.         break;
  922.     case eBlastx:
  923.     case eTblastx:
  924.         retval = 6;
  925.         break;
  926.     default:
  927.         NCBI_THROW(CBlastException, eBadParameter, 
  928.                    "Cannot get number of frames for invalid program type");
  929.     }
  931.     return retval;
  932. }
  933. END_SCOPE(blast)
  934. END_NCBI_SCOPE
  936. /* @} */
  938. /*
  939. * ===========================================================================
  940. *
  941. * $Log: blast_setup_cxx.cpp,v $
  942. * Revision 1000.7  2004/06/01 18:05:58  gouriano
  943. * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.68
  944. *
  945. * Revision 1.68  2004/05/21 21:41:02  gorelenk
  946. * Added PCH ncbi_pch.hpp
  947. *
  948. * Revision 1.67  2004/04/19 19:52:02  papadopo
  949. * correct off-by-one error in sequence size computation
  950. *
  951. * Revision 1.66  2004/04/16 14:28:49  papadopo
  952. * add use of eRPSBlast and eRPSTblastn programs
  953. *
  954. * Revision 1.65  2004/04/07 03:06:15  camacho
  955. * Added blast_encoding.[hc], refactoring blast_stat.[hc]
  956. *
  957. * Revision 1.64  2004/04/06 20:45:28  dondosha
  958. * Added function FindBlastDbPath: should be moved to seqdb library
  959. *
  960. * Revision 1.63  2004/03/24 19:14:48  dondosha
  961. * Fixed memory leaks
  962. *
  963. * Revision 1.62  2004/03/19 19:22:55  camacho
  964. * Move to doxygen group AlgoBlast, add missing CVS logs at EOF
  965. *
  966. * Revision 1.61  2004/03/15 19:57:52  dondosha
  967. * SetupSubjects takes just program argument instead of CBlastOptions*
  968. *
  969. * Revision 1.60  2004/03/11 17:26:46  dicuccio
  970. * Use NStr::ToUpper() instead of transform
  971. *
  972. * Revision 1.59  2004/03/09 18:53:25  dondosha
  973. * Do not set db length and number of sequences options to real values - these are calculated and assigned to parameters structure fields
  974. *
  975. * Revision 1.58  2004/03/06 00:39:47  camacho
  976. * Some refactorings, changed boolen parameter to enum in GetSequence
  977. *
  978. * Revision 1.57  2004/02/24 18:14:56  dondosha
  979. * Set the maximal length in the set of queries, when filling BlastQueryInfo structure
  980. *
  981. * Revision 1.56  2004/02/18 15:16:28  camacho
  982. * Consolidated ncbi2na mask definition
  983. *
  984. * Revision 1.55  2004/01/07 17:39:27  vasilche
  985. * Fixed include path to genbank loader.
  986. *
  987. * Revision 1.54  2003/12/29 17:00:57  camacho
  988. * Update comment
  989. *
  990. * Revision 1.53  2003/12/15 19:55:14  camacho
  991. * Minor fix to ensure exception safety
  992. *
  993. * Revision 1.52  2003/12/03 16:43:47  dondosha
  994. * Renamed BlastMask to BlastMaskLoc, BlastResults to BlastHSPResults
  995. *
  996. * Revision 1.51  2003/11/26 18:36:45  camacho
  997. * Renaming blast_option*pp -> blast_options*pp
  998. *
  999. * Revision 1.50  2003/11/26 18:24:00  camacho
  1000. * +Blast Option Handle classes
  1001. *
  1002. * Revision 1.49  2003/11/24 17:12:44  dondosha
  1003. * Query info structure does not have a total_length member any more; use last context offset
  1004. *
  1005. * Revision 1.48  2003/11/06 21:25:37  camacho
  1006. * Revert previous change, add assertions
  1007. *
  1008. * Revision 1.47  2003/11/04 17:14:22  dondosha
  1009. * Length in subject sequence block C structure should not include sentinels
  1010. *
  1011. * Revision 1.46  2003/10/31 19:45:03  camacho
  1012. * Fix setting of subject sequence length
  1013. *
  1014. * Revision 1.45  2003/10/31 16:53:32  dondosha
  1015. * Set length in query sequence block correctly
  1016. *
  1017. * Revision 1.44  2003/10/29 04:46:16  camacho
  1018. * Use fixed AutoPtr for GetSequence return value
  1019. *
  1020. * Revision 1.43  2003/10/27 21:27:36  camacho
  1021. * Remove extra argument to GetSequenceView, minor refactorings
  1022. *
  1023. * Revision 1.42  2003/10/22 14:21:55  camacho
  1024. * Added sanity checking assertions
  1025. *
  1026. * Revision 1.41  2003/10/21 13:04:54  camacho
  1027. * Fix bug in SetupSubjects, use sequence blk set functions
  1028. *
  1029. * Revision 1.40  2003/10/16 13:38:54  coulouri
  1030. * use anonymous exceptions to fix unreferenced variable compiler warning
  1031. *
  1032. * Revision 1.39  2003/10/15 18:18:00  camacho
  1033. * Fix to setup query info structure for proteins
  1034. *
  1035. * Revision 1.38  2003/10/15 15:09:32  camacho
  1036. * Changes from Mike DiCuccio to use GetSequenceView to retrieve sequences.
  1037. *
  1038. * Revision 1.37  2003/10/08 15:13:56  dondosha
  1039. * Test if subject mask is not NULL before converting to a C structure
  1040. *
  1041. * Revision 1.36  2003/10/08 15:05:47  dondosha
  1042. * Test if mask is not NULL before converting to a C structure
  1043. *
  1044. * Revision 1.35  2003/10/07 17:34:05  dondosha
  1045. * Add lower case masks to SSeqLocs forming the vector of sequence locations
  1046. *
  1047. * Revision 1.34  2003/10/03 16:12:18  dondosha
  1048. * Fix in previous change for plus strand search
  1049. *
  1050. * Revision 1.33  2003/10/02 22:10:46  dondosha
  1051. * Corrections for one-strand translated searches
  1052. *
  1053. * Revision 1.32  2003/09/30 03:23:18  camacho
  1054. * Fixes to FindMatrixPath
  1055. *
  1056. * Revision 1.31  2003/09/29 21:38:29  camacho
  1057. * Assign retval only when successfully found path to matrix
  1058. *
  1059. * Revision 1.30  2003/09/29 20:35:03  camacho
  1060. * Replace abort() with exception in GetNumberOfFrames
  1061. *
  1062. * Revision 1.29  2003/09/16 16:48:13  dondosha
  1063. * Use BLAST_PackDNA and BlastSetUp_SeqBlkNew from the core blast library for setting up subject sequences
  1064. *
  1065. * Revision 1.28  2003/09/12 17:52:42  camacho
  1066. * Stop using pair<> as return value from GetSequence
  1067. *
  1068. * Revision 1.27  2003/09/11 20:55:01  camacho
  1069. * Temporary fix for AutoPtr return value
  1070. *
  1071. * Revision 1.26  2003/09/11 17:45:03  camacho
  1072. * Changed CBlastOption -> CBlastOptions
  1073. *
  1074. * Revision 1.25  2003/09/10 04:27:43  camacho
  1075. * 1) Minor change to return type of GetSequence
  1076. * 2) Fix to previous revision
  1077. *
  1078. * Revision 1.24  2003/09/09 22:15:02  dondosha
  1079. * Added cast in return statement in GetSequence method, fixing compiler error
  1080. *
  1081. * Revision 1.23  2003/09/09 15:57:23  camacho
  1082. * Fix indentation
  1083. *
  1084. * Revision 1.22  2003/09/09 14:21:39  coulouri
  1085. * change blastkar.h to blast_stat.h
  1086. *
  1087. * Revision 1.21  2003/09/09 12:57:15  camacho
  1088. * + internal setup functions, use smart pointers to handle memory mgmt
  1089. *
  1090. * Revision 1.20  2003/09/05 19:06:31  camacho
  1091. * Use regular new to allocate genetic code string
  1092. *
  1093. * Revision 1.19  2003/09/03 19:36:27  camacho
  1094. * Fix include path for blast_setup.hpp
  1095. *
  1096. * Revision 1.18  2003/08/28 22:42:54  camacho
  1097. * Change BLASTGetSequence signature
  1098. *
  1099. * Revision 1.17  2003/08/28 15:49:02  madden
  1100. * Fix for packing DNA as well as correct buflen
  1101. *
  1102. * Revision 1.16  2003/08/25 16:24:14  camacho
  1103. * Updated BLASTGetMatrixPath
  1104. *
  1105. * Revision 1.15  2003/08/19 17:39:07  camacho
  1106. * Minor fix to use of metaregistry class
  1107. *
  1108. * Revision 1.14  2003/08/18 20:58:57  camacho
  1109. * Added blast namespace, removed *__.hpp includes
  1110. *
  1111. * Revision 1.13  2003/08/14 13:51:24  camacho
  1112. * Use CMetaRegistry class to load the ncbi config file
  1113. *
  1114. * Revision 1.12  2003/08/11 15:17:39  dondosha
  1115. * Added algo/blast/core to all #included headers
  1116. *
  1117. * Revision 1.11  2003/08/11 14:00:41  dicuccio
  1118. * Indenting changes.  Fixed use of C++ namespaces (USING_SCOPE(objects) inside of
  1119. * BEGIN_NCBI_SCOPE block)
  1120. *
  1121. * Revision 1.10  2003/08/08 19:43:07  dicuccio
  1122. * Compilation fixes: #include file rearrangement; fixed use of 'list' and
  1123. * 'vector' as variable names; fixed missing ostrea<< for __int64
  1124. *
  1125. * Revision 1.9  2003/08/04 15:18:23  camacho
  1126. * Minor fixes
  1127. *
  1128. * Revision 1.8  2003/08/01 22:35:02  camacho
  1129. * Added function to get matrix path (fixme)
  1130. *
  1131. * Revision 1.7  2003/08/01 17:40:56  dondosha
  1132. * Use renamed functions and structures from local blastkar.h
  1133. *
  1134. * Revision 1.6  2003/07/31 19:45:33  camacho
  1135. * Eliminate Ptr notation
  1136. *
  1137. * Revision 1.5  2003/07/30 15:00:01  camacho
  1138. * Do not use Malloc/MemNew/MemFree
  1139. *
  1140. * Revision 1.4  2003/07/25 13:55:58  camacho
  1141. * Removed unnecessary #includes
  1142. *
  1143. * Revision 1.3  2003/07/24 18:22:50  camacho
  1144. * #include blastkar.h
  1145. *
  1146. * Revision 1.2  2003/07/23 21:29:06  camacho
  1147. * Update BLASTFindGeneticCode to get genetic code string with C++ toolkit
  1148. *
  1149. * Revision 1.1  2003/07/10 18:34:19  camacho
  1150. * Initial revision
  1151. *
  1152. *
  1153. * ===========================================================================
  1154. */