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

生物技术

开发平台:

C/C++

  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: struct_dp_demo.cpp,v $
  4.  * PRODUCTION Revision 1000.2  2004/06/01 18:11:14  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.11
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  9. /*  $Id: struct_dp_demo.cpp,v 1000.2 2004/06/01 18:11:14 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:  Paul Thiessen
  35. *
  36. * File Description:
  37. *      Test/demo for dynamic programming-based alignment algorithms
  38. *
  39. * ===========================================================================
  40. */
  41. #include <ncbi_pch.hpp>
  42. #include <corelib/ncbistd.hpp>
  43. #include <corelib/ncbi_limits.h>
  44. #include <util/tables/raw_scoremat.h>
  45. #include <map>
  46. #include <string>
  47. #include <algo/structure/struct_dp/struct_dp.h>
  48. #include "struct_dp_demo.hpp"
  49. USING_NCBI_SCOPE;
  50. BEGIN_SCOPE(struct_dp)
  51. #define ERROR_MESSAGE(s) ERR_POST(Error << "struct_dp_demo: " << s << '!')
  52. #define INFO_MESSAGE(s) ERR_POST(Info << "struct_dp_demo: " << s)
  53. void DPApp::Init(void)
  54. {
  55. }
  56. static inline char ScreenResidueCharacter(char original)
  57. {
  58.     char ch = toupper(original);
  59.     switch (ch) {
  60.         case 'A': case 'R': case 'N': case 'D': case 'C':
  61.         case 'Q': case 'E': case 'G': case 'H': case 'I':
  62.         case 'L': case 'K': case 'M': case 'F': case 'P':
  63.         case 'S': case 'T': case 'W': case 'Y': case 'V':
  64.         case 'B': case 'Z':
  65.             break;
  66.         default:
  67.             ch = 'X'; // make all but natural aa's just 'X'
  68.     }
  69.     return ch;
  70. }
  71. static int GetBLOSUM62Score(char a, char b)
  72. {
  73.     static SNCBIFullScoreMatrix Blosum62Matrix;
  74.     static bool unpacked = false;
  75.     if (!unpacked) {
  76.         NCBISM_Unpack(&NCBISM_Blosum62, &Blosum62Matrix);
  77.         unpacked = true;
  78.     }
  79.     return Blosum62Matrix.s[ScreenResidueCharacter(a)][ScreenResidueCharacter(b)];
  80. }
  81. static DP_BlockInfo *blocks = NULL;
  82. static string subject, query;
  83. extern "C" {
  84. int ScoreByBlosum62(unsigned int block, unsigned int queryPos)
  85. {
  86.     if (!blocks || block >= blocks->nBlocks ||
  87.         queryPos + blocks->blockSizes[block] > query.size())
  88.     {
  89.         ERROR_MESSAGE("ScoreByBlosum62() - invalid parameters");
  90.         return DP_NEGATIVE_INFINITY;
  91.     }
  92.     int score = 0;
  93.     for (int i=0; i<blocks->blockSizes[block]; i++)
  94.         score += GetBLOSUM62Score(
  95.             subject[blocks->blockPositions[block] + i],
  96.             query[queryPos + i]);
  97. //    INFO_MESSAGE("score for block " << block << " pos " << queryPos << ": " << score);
  98.     return score;
  99. }
  100. } // extern "C"
  101. int DPApp::Run(void)
  102. {
  103.     // this demo tries to reproduce the VAST alignment of 1DOI and 1FRD
  104.     subject = "PTVEYLNYEVVDDNGWDMYDDDVFGEASDMDLDDEDYGSLEVNEGEYILEAAEAQGYDWPFSC"
  105.         "RAGACANCAAIVLEGDIDMDMQQILSDEEVEDKNVRLTCIGSPDADEVKIVYNAKHLDYLQNRVI";
  106.     blocks = DP_CreateBlockInfo(5);
  107.     blocks->blockPositions[0] = 0;
  108.     blocks->blockSizes[0] = 7;
  109.     blocks->maxLoops[0] = 100;
  110.     blocks->blockPositions[1] = 35;
  111.     blocks->blockSizes[1] = 49;
  112.     blocks->maxLoops[1] = 100;
  113.     blocks->blockPositions[2] = 86;
  114.     blocks->blockSizes[2] = 8;
  115.     blocks->maxLoops[2] = 100;
  116.     blocks->blockPositions[3] = 97;
  117.     blocks->blockSizes[3] = 10;
  118.     blocks->maxLoops[3] = 100;
  119. //    blocks->maxLoops[3] = 0;
  120. //    blocks->freezeBlocks[3] = 74;
  121.     blocks->blockPositions[4] = 109;
  122.     blocks->blockSizes[4] = 11;
  123. //    blocks->freezeBlocks[4] = 85;
  124. //    for (int b=0; b<blocks->nBlocks; b++)
  125. //        INFO_MESSAGE("Block " << (b+1) << ": "
  126. //            << subject.substr(blocks->blockPositions[b], blocks->blockSizes[b]));
  127.     query = "ASYQVRLINKKQDIDTTIEIDEETTILDGAEENGIELPFSCHSGSCSSCVGKVVEGEVDQSDQ"
  128.                 "IFLDDEQMGKGFALLCVTYPRSNCTIKTHQEPYLA";
  129.     DP_AlignmentResult *alignment;
  130.     int result = DP_LocalBlockAlign(blocks, ScoreByBlosum62, 0, query.size() - 1, &alignment);
  131.     // crudely print out alignment result
  132.     if (result == STRUCT_DP_FOUND_ALIGNMENT) {
  133.         INFO_MESSAGE("Found alignment, total score: " << alignment->score
  134.             << " for " << alignment->nBlocks << " blocks");
  135.         unsigned int block;
  136.         for (block=0; block<alignment->nBlocks; block++) {
  137.             unsigned int
  138.                 subjectStart = blocks->blockPositions[block + alignment->firstBlock],
  139.                 queryStart = alignment->blockPositions[block],
  140.                 blockSize = blocks->blockSizes[block + alignment->firstBlock];
  141.             INFO_MESSAGE(
  142.                 "Block " << (block + alignment->firstBlock + 1) << ", score "
  143.                 << ScoreByBlosum62(block + alignment->firstBlock, queryStart) << ":n"
  144.                 << "S: " << subject.substr(subjectStart, blockSize)
  145.                 << ' ' << (subjectStart + 1) << '-' << (subjectStart + blockSize) << 'n'
  146.                 << "Q: " << query.substr(queryStart, blockSize)
  147.                 << ' ' << (queryStart + 1) << '-' << (queryStart + blockSize)
  148.             );
  149.         }
  150.         DP_DestroyAlignmentResult(alignment);
  151.     }
  152.     else if (result == STRUCT_DP_NO_ALIGNMENT) {
  153.         INFO_MESSAGE("Found no significant alignment");
  154.     }
  155.     else if (result == STRUCT_DP_PARAMETER_ERROR || result == STRUCT_DP_ALGORITHM_ERROR) {
  156.         ERROR_MESSAGE("DP_GlobalBlockAlign() failed");
  157.     }
  158.     DP_DestroyBlockInfo(blocks);
  159.     return 0;
  160. }
  161. END_SCOPE(struct_dp)
  162. USING_NCBI_SCOPE;
  163. USING_SCOPE(struct_dp);
  164. int main(int argc, const char* argv[])
  165. {
  166.     SetDiagStream(&NcbiCerr); // send all diagnostic messages to cerr
  167.     SetDiagPostLevel(eDiag_Info);   // show all messages
  168.     DPApp app;
  169.     return app.AppMain(argc, argv, NULL, eDS_Default, NULL);    // don't use config file
  170. }
  171. /*
  172. * ---------------------------------------------------------------------------
  173. * $Log: struct_dp_demo.cpp,v $
  174. * Revision 1000.2  2004/06/01 18:11:14  gouriano
  175. * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.11
  176. *
  177. * Revision 1.11  2004/05/21 21:41:04  gorelenk
  178. * Added PCH ncbi_pch.hpp
  179. *
  180. * Revision 1.10  2004/02/19 02:21:20  thiessen
  181. * fix struct_dp paths
  182. *
  183. * Revision 1.9  2003/12/11 19:47:48  thiessen
  184. * use built-in blosum62
  185. *
  186. * Revision 1.8  2003/07/11 15:27:48  thiessen
  187. * add DP_ prefix to globals
  188. *
  189. * Revision 1.7  2003/06/22 12:11:38  thiessen
  190. * add local alignment
  191. *
  192. * Revision 1.6  2003/06/18 22:25:37  thiessen
  193. * simplify alignment printout code, yet again... ;)
  194. *
  195. * Revision 1.5  2003/06/18 22:11:00  thiessen
  196. * simplify alignment printout code, again... ;)
  197. *
  198. * Revision 1.4  2003/06/18 22:06:03  thiessen
  199. * simplify alignment printout code
  200. *
  201. * Revision 1.3  2003/06/18 21:46:10  thiessen
  202. * add traceback, alignment return structure
  203. *
  204. * Revision 1.2  2003/06/18 19:10:17  thiessen
  205. * fix lf issues
  206. *
  207. * Revision 1.1  2003/06/18 16:54:29  thiessen
  208. * initial checkin; working global block aligner and demo app
  209. *
  210. */