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

生物技术

开发平台:

C/C++

  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: blast_input.cpp,v $
  4.  * PRODUCTION Revision 1000.3  2004/06/01 18:06:38  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.22
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  9. /*  $Id: blast_input.cpp,v 1000.3 2004/06/01 18:06:38 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 official 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.  * Authors:  Christiam Camacho
  35.  *
  36.  */
  37. /** @file blast_input.cpp
  38.  * Reading FASTA from an input file
  39.  */
  40. static char const rcsid[] = 
  41.     "$Id: blast_input.cpp,v 1000.3 2004/06/01 18:06:38 gouriano Exp $";
  42. #include <ncbi_pch.hpp>
  43. #include <serial/iterator.hpp>
  44. #include <objmgr/util/sequence.hpp>
  45. #include <objtools/readers/fasta.hpp>
  46. #include <objtools/readers/reader_exception.hpp>
  47. #include <objects/seq/Bioseq.hpp>
  48. #include <objects/seqloc/Seq_loc.hpp>
  49. #include <objects/seqloc/Seq_interval.hpp>
  50. #include <objmgr/scope.hpp>
  51. #include <algo/blast/api/blast_aux.hpp>
  52. #include "blast_input.hpp"
  53. BEGIN_NCBI_SCOPE
  54. USING_SCOPE(objects);
  55. BEGIN_SCOPE(blast)
  56. TSeqLocVector
  57. BLASTGetSeqLocFromStream(CNcbiIstream& in, CObjectManager& objmgr, 
  58.                          ENa_strand strand, TSeqPos from, TSeqPos to, 
  59.                          int *counter, bool get_lcase_mask)
  60. {
  61.     TSeqLocVector retval;
  62.     CRef<CSeq_entry> seq_entry;
  63.     vector<CConstRef<CSeq_loc> > lcase_mask;
  64.     CRef<CScope> scope(new CScope(objmgr));
  65.     scope->AddDefaults();
  66.     if (get_lcase_mask) {
  67.         if ( !(seq_entry = ReadFasta(in, fReadFasta_AllSeqIds, counter, 
  68.                                      &lcase_mask)))
  69.             NCBI_THROW(CObjReaderException, eInvalid, 
  70.                        "Could not retrieve seq entry");
  71.     } else {
  72.         if ( !(seq_entry = ReadFasta(in, fReadFasta_AllSeqIds, counter)))
  73.             NCBI_THROW(CObjReaderException, eInvalid, 
  74.                        "Could not retrieve seq entry");
  75.     }
  76.     int index = 0;
  77.     scope->AddTopLevelSeqEntry(*seq_entry);
  78.     for (CTypeConstIterator<CBioseq> itr(ConstBegin(*seq_entry)); itr; ++itr) {
  79.         CRef<CSeq_loc> seqloc(new CSeq_loc());
  80.         TSeqPos seq_length = sequence::GetLength(*itr->GetId().front(), 
  81.                                                 scope)-1;
  82.         if (to > 0 && to < seq_length)
  83.             seqloc->SetInt().SetTo(to);
  84.         else
  85.             seqloc->SetInt().SetTo(seq_length);
  86.         if (from > 0 && from < seq_length && from < to)
  87.             seqloc->SetInt().SetFrom(from);
  88.         else
  89.             seqloc->SetInt().SetFrom(0);
  90.         seqloc->SetInt().SetStrand(strand);
  91.         seqloc->SetInt().SetId().Assign(*itr->GetId().front());
  92.         //CRef<CScope> s(scope);
  93.         SSeqLoc sl(seqloc, scope);
  94.         if (get_lcase_mask) {
  95.             sl.mask.Reset(lcase_mask[index++]);
  96.         }
  97.         retval.push_back(sl);
  98.     }
  99.     return retval;
  100. }
  101. END_SCOPE(blast)
  102. END_NCBI_SCOPE