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

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

生物技术

开发平台：

C/C++

cn3d_blast.cpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: cn3d_blast.cpp,v $
  4.  * PRODUCTION Revision 1000.2  2004/06/01 18:28:09  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.34
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /*  $Id: cn3d_blast.cpp,v 1000.2 2004/06/01 18:28:09 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. * Authors:  Paul Thiessen
  36. *
  37. * File Description:
  38. *      module for aligning with BLAST and related algorithms
  39. *
  40. * ===========================================================================
  41. */
  43. #ifdef _MSC_VER
  44. #pragma warning(disable:4018)   // disable signed/unsigned mismatch warning in MSVC
  45. #endif
  47. #include <ncbi_pch.hpp>
  48. #include <corelib/ncbistd.hpp>
  50. #include <objects/seqalign/Seq_align.hpp>
  51. #include <objects/seqalign/Dense_seg.hpp>
  52. #include <objects/seqalign/Score.hpp>
  53. #include <objects/general/Object_id.hpp>
  55. #ifdef __WXMSW__
  56. #include <windows.h>
  57. #include <wx/msw/winundef.h>
  58. #endif
  59. #include <wx/wx.h>
  61. // C stuff
  62. #include <objseq.h>
  63. #include <blast.h>
  64. #include <blastkar.h>
  66. #include "cn3d_blast.hpp"
  67. #include "structure_set.hpp"
  68. #include "sequence_set.hpp"
  69. #include "block_multiple_alignment.hpp"
  70. #include "cn3d_tools.hpp"
  71. #include "asn_converter.hpp"
  72. #include "molecule_identifier.hpp"
  73. #include "cn3d_threader.hpp"
  75. // hack so I can catch memory leaks specific to this module, at the line where allocation occurs
  76. #ifdef _DEBUG
  77. #ifdef MemNew
  78. #undef MemNew
  79. #endif
  80. #define MemNew(sz) memset(malloc(sz), 0, (sz))
  81. #endif
  83. USING_NCBI_SCOPE;
  84. USING_SCOPE(objects);
  87. BEGIN_SCOPE(Cn3D)
  89. const string BLASTer::BLASTResidues = "-ABCDEFGHIKLMNPQRSTVWXYZU*";
  91. // gives BLAST residue number for a character (or # for 'X' if char not found)
  92. int LookupBLASTResidueNumberFromCharacter(char r)
  93. {
  94.     typedef map < char, int > Char2Int;
  95.     static Char2Int charMap;
  97.     if (charMap.size() == 0) {
  98.         for (int i=0; i<BLASTer::BLASTResidues.size(); ++i)
  99.             charMap[BLASTer::BLASTResidues[i]] = i;
  100.     }
  102.     Char2Int::const_iterator n = charMap.find(toupper(r));
  103.     if (n != charMap.end())
  104.         return n->second;
  105.     else
  106.         return charMap.find('X')->second;
  107. }
  109. // gives BLAST residue number for a threader residue number (or # for 'X' if char == -1)
  110. int LookupBLASTResidueNumberFromThreaderResidueNumber(char r)
  111. {
  112.     r = toupper(r);
  113.     return LookupBLASTResidueNumberFromCharacter(
  114.             (r >= 0 && r < Threader::ThreaderResidues.size()) ? Threader::ThreaderResidues[r] : 'X');
  115. }
  117. #ifdef _DEBUG
  118. #define PRINT_PSSM // for testing/debugging
  119. #endif
  121. static BLAST_OptionsBlkPtr CreateBlastOptionsBlk(void)
  122. {
  123.     BLAST_OptionsBlkPtr bob = BLASTOptionNew("blastp", true);
  124.     bob->db_length = 2717223;    // size of CDD database v1.62
  125.     bob->scalingFactor = 1.0;
  126.     return bob;
  127. }
  129. static void SetEffectiveSearchSpaceSize(BLAST_OptionsBlkPtr bob, Int4 queryLength)
  130. {
  131.     bob->searchsp_eff = BLASTCalculateSearchSpace(bob, 1, bob->db_length, queryLength);
  132. }
  134. void BLASTer::CreateNewPairwiseAlignmentsByBlast(const BlockMultipleAlignment *multiple,
  135.     const AlignmentList& toRealign, AlignmentList *newAlignments, bool usePSSM)
  136. {
  137.     if (usePSSM && !multiple) {
  138.         ERRORMSG("usePSSM true, but NULL multiple alignment");
  139.         return;
  140.     }
  142.     // set up BLAST stuff
  143.     BLAST_OptionsBlkPtr options = CreateBlastOptionsBlk();
  144.     if (!options) {
  145.         ERRORMSG("BLASTOptionNew() failed");
  146.         return;
  147.     }
  149.     // create Seq-loc intervals
  150.     SeqLocPtr masterSeqLoc = (SeqLocPtr) MemNew(sizeof(SeqLoc));
  151.     masterSeqLoc->choice = SEQLOC_INT;
  152.     SeqIntPtr masterSeqInt = SeqIntNew();
  153.     masterSeqLoc->data.ptrvalue = masterSeqInt;
  154.     if (multiple) {
  155.         BlockMultipleAlignment::UngappedAlignedBlockList uaBlocks;
  156.         multiple->GetUngappedAlignedBlocks(&uaBlocks);
  157.         if (uaBlocks.size() > 0) {
  158.             int excess = 0;
  159.             if (!RegistryGetInteger(REG_ADVANCED_SECTION, REG_FOOTPRINT_RES, &excess))
  160.                 WARNINGMSG("Can't get footprint excess residues from registry");
  161.             masterSeqInt->from = uaBlocks.front()->GetRangeOfRow(0)->from - excess;
  162.             if (masterSeqInt->from < 0)
  163.                 masterSeqInt->from = 0;
  164.             masterSeqInt->to = uaBlocks.back()->GetRangeOfRow(0)->to + excess;
  165.             if (masterSeqInt->to >= multiple->GetMaster()->Length())
  166.                 masterSeqInt->to = multiple->GetMaster()->Length() - 1;
  167.         } else {
  168.             masterSeqInt->from = 0;
  169.             masterSeqInt->to = multiple->GetMaster()->Length() - 1;
  170.         }
  171.     }
  172.     SeqLocPtr slaveSeqLoc = (SeqLocPtr) MemNew(sizeof(SeqLoc));
  173.     slaveSeqLoc->choice = SEQLOC_INT;
  174.     SeqIntPtr slaveSeqInt = SeqIntNew();
  175.     slaveSeqLoc->data.ptrvalue = slaveSeqInt;
  177.     // create BLAST_Matrix if using PSSM from multiple
  178.     const BLAST_Matrix *BLASTmatrix = NULL;
  179.     if (usePSSM) BLASTmatrix = multiple->GetPSSM();
  181.     string err;
  182.     AlignmentList::const_iterator s, se = toRealign.end();
  183.     for (s=toRealign.begin(); s!=se; ++s) {
  184.         if (multiple && (*s)->GetMaster() != multiple->GetMaster())
  185.             ERRORMSG("master sequence mismatch");
  187.         // get C Bioseq for master
  188.         const Sequence *masterSeq = multiple ? multiple->GetMaster() : (*s)->GetMaster();
  189.         Bioseq *masterBioseq = masterSeq->parentSet->GetOrCreateBioseq(masterSeq);
  190.         masterSeqInt->id = masterBioseq->id;
  192.         // override master alignment interval if specified
  193.         if (!multiple) {
  194.             if ((*s)->alignMasterFrom >= 0 && (*s)->alignMasterFrom < masterSeq->Length() &&
  195.                 (*s)->alignMasterTo >= 0 && (*s)->alignMasterTo < masterSeq->Length() &&
  196.                 (*s)->alignMasterFrom <= (*s)->alignMasterTo)
  197.             {
  198.                 masterSeqInt->from = (*s)->alignMasterFrom;
  199.                 masterSeqInt->to = (*s)->alignMasterTo;
  200.             } else {
  201.                 masterSeqInt->from = 0;
  202.                 masterSeqInt->to = masterSeq->Length() - 1;
  203.             }
  204.         }
  206.         // get C Bioseq for slave of each incoming (pairwise) alignment
  207.         const Sequence *slaveSeq = (*s)->GetSequenceOfRow(1);
  208.         Bioseq *slaveBioseq = slaveSeq->parentSet->GetOrCreateBioseq(slaveSeq);
  210.         // setup Seq-loc interval for slave
  211.         slaveSeqInt->id = slaveBioseq->id;
  212.         slaveSeqInt->from =
  213.             ((*s)->alignSlaveFrom >= 0 && (*s)->alignSlaveFrom < slaveSeq->Length()) ?
  214.                 (*s)->alignSlaveFrom : 0;
  215.         slaveSeqInt->to =
  216.             ((*s)->alignSlaveTo >= 0 && (*s)->alignSlaveTo < slaveSeq->Length()) ?
  217.                 (*s)->alignSlaveTo : slaveSeq->Length() - 1;
  218.         INFOMSG("for BLAST: master range " <<
  219.                 (masterSeqInt->from + 1) << " to " << (masterSeqInt->to + 1) << ", slave range " <<
  220.                 (slaveSeqInt->from + 1) << " to " << (slaveSeqInt->to + 1));
  222.         // set search space size
  223.         SetEffectiveSearchSpaceSize(options, slaveSeqInt->to - slaveSeqInt->from + 1);
  224.         INFOMSG("effective search space size: " << ((int) options->searchsp_eff));
  226.         // actually do the BLAST alignment
  227.         SeqAlign *salp = NULL;
  228.         if (usePSSM) {
  229.             INFOMSG("calling BlastTwoSequencesByLocWithCallback(); PSSM db_length "
  230.                 << (int) options->db_length << ", K " << BLASTmatrix->karlinK);
  231.             salp = BlastTwoSequencesByLocWithCallback(
  232.                 masterSeqLoc, slaveSeqLoc,
  233.                 "blastp", options,
  234.                 NULL, NULL, NULL,
  235.                 const_cast<BLAST_Matrix*>(BLASTmatrix));
  236.         } else {
  237.             INFOMSG("calling BlastTwoSequencesByLoc()");
  238.             salp = BlastTwoSequencesByLoc(masterSeqLoc, slaveSeqLoc, "blastp", options);
  239.         }
  241.         // create new alignment structure
  242.         BlockMultipleAlignment::SequenceList *seqs = new BlockMultipleAlignment::SequenceList(2);
  243.         (*seqs)[0] = masterSeq;
  244.         (*seqs)[1] = slaveSeq;
  245.         BlockMultipleAlignment *newAlignment =
  246.             new BlockMultipleAlignment(seqs, slaveSeq->parentSet->alignmentManager);
  247.         newAlignment->SetRowDouble(0, kMax_Double);
  248.         newAlignment->SetRowDouble(1, kMax_Double);
  250.         // process the result
  251.         if (!salp) {
  252.             WARNINGMSG("BLAST failed to find a significant alignment");
  253.         } else {
  255.             // convert C SeqAlign to C++ for convenience
  256. #ifdef _DEBUG
  257.             AsnIoPtr aip = AsnIoOpen("seqalign.txt", "w");
  258.             SeqAlignAsnWrite(salp, aip, NULL);
  259.             AsnIoFree(aip, true);
  260. #endif
  261.             CSeq_align sa;
  262.             bool okay = ConvertAsnFromCToCPP(salp, (AsnWriteFunc) SeqAlignAsnWrite, &sa, &err);
  263.             SeqAlignSetFree(salp);
  264.             if (!okay) {
  265.                 ERRORMSG("Conversion of SeqAlign to C++ object failed");
  266.                 continue;
  267.             }
  269.             // make sure the structure of this SeqAlign is what we expect
  270.             if (!sa.IsSetDim() || sa.GetDim() != 2 ||
  271.                 !sa.GetSegs().IsDenseg() || sa.GetSegs().GetDenseg().GetDim() != 2 ||
  272.                 !masterSeq->identifier->MatchesSeqId(*(sa.GetSegs().GetDenseg().GetIds().front())) ||
  273.                 !slaveSeq->identifier->MatchesSeqId(*(sa.GetSegs().GetDenseg().GetIds().back()))) {
  274.                 ERRORMSG("Confused by BLAST result format");
  275.                 continue;
  276.             }
  278.             // fill out with blocks from BLAST alignment
  279.             CDense_seg::TStarts::const_iterator iStart = sa.GetSegs().GetDenseg().GetStarts().begin();
  280.             CDense_seg::TLens::const_iterator iLen = sa.GetSegs().GetDenseg().GetLens().begin();
  281. #ifdef _DEBUG
  282.             int raw_pssm = 0, raw_bl62 = 0;
  283. #endif
  284.             for (int i=0; i<sa.GetSegs().GetDenseg().GetNumseg(); ++i) {
  285.                 int masterStart = *(iStart++), slaveStart = *(iStart++), len = *(iLen++);
  286.                 if (masterStart >= 0 && slaveStart >= 0 && len > 0) {
  287.                     UngappedAlignedBlock *newBlock = new UngappedAlignedBlock(newAlignment);
  288.                     newBlock->SetRangeOfRow(0, masterStart, masterStart + len - 1);
  289.                     newBlock->SetRangeOfRow(1, slaveStart, slaveStart + len - 1);
  290.                     newBlock->width = len;
  291.                     newAlignment->AddAlignedBlockAtEnd(newBlock);
  293. #ifdef _DEBUG
  294.                     // calculate score manually, to compare with that returned by BLAST
  295.                     for (int j=0; j<len; ++j) {
  296.                         if (usePSSM)
  297.                             raw_pssm += BLASTmatrix->matrix
  298.                                 [masterStart + j]
  299.                                 [LookupBLASTResidueNumberFromCharacter(
  300.                                     slaveSeq->sequenceString[slaveStart + j])];
  301.                         raw_bl62 += GetBLOSUM62Score(
  302.                             masterSeq->sequenceString[masterStart + j],
  303.                             slaveSeq->sequenceString[slaveStart + j]);
  304.                     }
  305. #endif
  306.                 }
  307. #ifdef _DEBUG
  308.                 else if ((masterStart < 0 || slaveStart < 0) && len > 0) {
  309.                     int gap = options->gap_open + options->gap_extend * len;
  310.                     if (usePSSM) raw_pssm -= gap;
  311.                     raw_bl62 -= gap;
  312.                 }
  313. #endif
  314.             }
  315. #ifdef _DEBUG
  316.             TRACEMSG("Calculated raw score by PSSM: " << raw_pssm);
  317.             TRACEMSG("Calculated raw score by BLOSUM62: " << raw_bl62);
  318. #endif
  320.             // get scores
  321.             if (sa.IsSetScore()) {
  322.                 wxString scores;
  323.                 scores.Printf("BLAST result for %s vs. %s:",
  324.                     masterSeq->identifier->ToString().c_str(),
  325.                     slaveSeq->identifier->ToString().c_str());
  327.                 CSeq_align::TScore::const_iterator sc, sce = sa.GetScore().end();
  328.                 for (sc=sa.GetScore().begin(); sc!=sce; ++sc) {
  329.                     if ((*sc)->IsSetId() && (*sc)->GetId().IsStr()) {
  331.                         // E-value (put in status line and double values)
  332.                         if ((*sc)->GetValue().IsReal() && (*sc)->GetId().GetStr() == "e_value") {
  333.                             newAlignment->SetRowDouble(0, (*sc)->GetValue().GetReal());
  334.                             newAlignment->SetRowDouble(1, (*sc)->GetValue().GetReal());
  335.                             wxString status;
  336.                             status.Printf("E-value %g", (*sc)->GetValue().GetReal());
  337.                             newAlignment->SetRowStatusLine(0, status.c_str());
  338.                             newAlignment->SetRowStatusLine(1, status.c_str());
  339.                             wxString tmp = scores;
  340.                             scores.Printf("%s E-value: %g", tmp.c_str(), (*sc)->GetValue().GetReal());
  341.                         }
  343.                         // raw score
  344.                         if ((*sc)->GetValue().IsInt() && (*sc)->GetId().GetStr() == "score") {
  345.                             wxString tmp = scores;
  346.                             scores.Printf("%s raw: %i", tmp.c_str(), (*sc)->GetValue().GetInt());
  347.                         }
  349.                         // bit score
  350.                         if ((*sc)->GetValue().IsReal() && (*sc)->GetId().GetStr() == "bit_score") {
  351.                             wxString tmp = scores;
  352.                             scores.Printf("%s bit score: %g", tmp.c_str(), (*sc)->GetValue().GetReal());
  353.                         }
  354.                     }
  355.                 }
  357.                 INFOMSG(scores.c_str());
  358.             }
  359.         }
  361.         // finalize and and add new alignment to list
  362.         if (newAlignment->AddUnalignedBlocks() && newAlignment->UpdateBlockMapAndColors(false)) {
  363.             newAlignments->push_back(newAlignment);
  364.         } else {
  365.             ERRORMSG("error finalizing alignment");
  366.             delete newAlignment;
  367.         }
  368.     }
  370.     BLASTOptionDelete(options);
  371.     masterSeqInt->id = NULL;    // don't free Seq-id, since it belongs to the Bioseq
  372.     SeqIntFree(masterSeqInt);
  373.     MemFree(masterSeqLoc);
  374.     slaveSeqInt->id = NULL;
  375.     SeqIntFree(slaveSeqInt);
  376.     MemFree(slaveSeqLoc);
  377. }
  379. void BLASTer::CalculateSelfHitScores(const BlockMultipleAlignment *multiple)
  380. {
  381.     if (!multiple) {
  382.         ERRORMSG("NULL multiple alignment");
  383.         return;
  384.     }
  386.     // set up BLAST stuff
  387.     BLAST_OptionsBlkPtr options = CreateBlastOptionsBlk();
  388.     if (!options) {
  389.         ERRORMSG("BLASTOptionNew() failed");
  390.         return;
  391.     }
  393.     // get master sequence
  394.     const Sequence *masterSeq = multiple->GetMaster();
  395.     Bioseq *masterBioseq = masterSeq->parentSet->GetOrCreateBioseq(masterSeq);
  397.     // create Seq-loc interval for master and slave
  398.     SeqLocPtr masterSeqLoc = (SeqLocPtr) MemNew(sizeof(SeqLoc));
  399.     masterSeqLoc->choice = SEQLOC_INT;
  400.     SeqIntPtr masterSeqInt = SeqIntNew();
  401.     masterSeqLoc->data.ptrvalue = masterSeqInt;
  402.     masterSeqInt->id = masterBioseq->id;
  403.     BlockMultipleAlignment::UngappedAlignedBlockList uaBlocks;
  404.     multiple->GetUngappedAlignedBlocks(&uaBlocks);
  405.     if (uaBlocks.size() == 0) {
  406.         ERRORMSG("Self-hit requires at least one aligned block");
  407.         return;
  408.     }
  409.     int excess = 0;
  410.     if (!RegistryGetInteger(REG_ADVANCED_SECTION, REG_FOOTPRINT_RES, &excess))
  411.         WARNINGMSG("Can't get footprint excess residues from registry");
  412.     masterSeqInt->from = uaBlocks.front()->GetRangeOfRow(0)->from - excess;
  413.     if (masterSeqInt->from < 0) masterSeqInt->from = 0;
  414.     masterSeqInt->to = uaBlocks.back()->GetRangeOfRow(0)->to + excess;
  415.     if (masterSeqInt->to >= masterSeq->Length()) masterSeqInt->to = masterSeq->Length() - 1;
  417.     SeqLocPtr slaveSeqLoc = (SeqLocPtr) MemNew(sizeof(SeqLoc));
  418.     slaveSeqLoc->choice = SEQLOC_INT;
  419.     SeqIntPtr slaveSeqInt = SeqIntNew();
  420.     slaveSeqLoc->data.ptrvalue = slaveSeqInt;
  422.     // create BLAST_Matrix if using PSSM from multiple
  423.     const BLAST_Matrix *BLASTmatrix = multiple->GetPSSM();
  425.     string err;
  426.     int row;
  427.     for (row=0; row<multiple->NRows(); ++row) {
  429.         // get C Bioseq for each slave
  430.         const Sequence *slaveSeq = multiple->GetSequenceOfRow(row);
  431.         Bioseq *slaveBioseq = slaveSeq->parentSet->GetOrCreateBioseq(slaveSeq);
  433.         // setup Seq-loc interval for slave
  434.         slaveSeqInt->id = slaveBioseq->id;
  435.         slaveSeqInt->from = uaBlocks.front()->GetRangeOfRow(row)->from - excess;
  436.         if (slaveSeqInt->from < 0) slaveSeqInt->from = 0;
  437.         slaveSeqInt->to = uaBlocks.back()->GetRangeOfRow(row)->to + excess;
  438.         if (slaveSeqInt->to >= slaveSeq->Length()) slaveSeqInt->to = slaveSeq->Length() - 1;
  439. //        TESTMSG("for BLAST: master range " <<
  440. //                (masterSeqInt->from + 1) << " to " << (masterSeqInt->to + 1) << ", slave range " <<
  441. //                (slaveSeqInt->from + 1) << " to " << (slaveSeqInt->to + 1));
  443.         // set search space size
  444.         SetEffectiveSearchSpaceSize(options, slaveSeqInt->to - slaveSeqInt->from + 1);
  446.         // actually do the BLAST alignment
  448. //        TRACEMSG("calling BlastTwoSequencesByLocWithCallback()");
  449.         SeqAlign *salp = BlastTwoSequencesByLocWithCallback(
  450.             masterSeqLoc, slaveSeqLoc,
  451.             "blastp", options,
  452.             NULL, NULL, NULL,
  453.             const_cast<BLAST_Matrix*>(BLASTmatrix));
  455. //        TRACEMSG("calling B2SPssmMultipleQueries()");
  456. //        SeqLocPtr target[1];
  457. //        target[0] = slaveSeqLoc;
  458. //        options->searchsp_eff = 0;
  459. //        SeqAlignPtr *psalp = B2SPssmMultipleQueries(masterSeqLoc,
  460. //            const_cast<BLAST_Matrix*>(BLASTmatrix), target, 1, options);
  461. //        SeqAlign *salp = psalp[0];
  463.         // process the result
  464.         double score = -1.0;
  465.         if (salp) {
  466.             // convert C SeqAlign to C++ for convenience
  467.             CSeq_align sa;
  468.             bool okay = ConvertAsnFromCToCPP(salp, (AsnWriteFunc) SeqAlignAsnWrite, &sa, &err);
  469.             SeqAlignSetFree(salp);
  470.             if (!okay) {
  471.                 ERRORMSG("Conversion of SeqAlign to C++ object failed");
  472.                 continue;
  473.             }
  475.             // get score
  476.             if (sa.IsSetScore()) {
  477.                 CSeq_align::TScore::const_iterator sc, sce = sa.GetScore().end();
  478.                 for (sc=sa.GetScore().begin(); sc!=sce; ++sc) {
  479.                     if ((*sc)->GetValue().IsReal() && (*sc)->IsSetId() &&
  480.                         (*sc)->GetId().IsStr() && (*sc)->GetId().GetStr() == "e_value") {
  481.                         score = (*sc)->GetValue().GetReal();
  482.                         break;
  483.                     }
  484.                 }
  485.             }
  486.             if (score < 0.0) ERRORMSG("Got back Seq-align with no E-value");
  488.         }
  490.         // set score in row
  491.         multiple->SetRowDouble(row, score);
  492.         wxString status;
  493.         if (score >= 0.0)
  494.             status.Printf("Self hit E-value: %g", score);
  495.         else
  496.             status = "No detectable self hit";
  497.         multiple->SetRowStatusLine(row, status.c_str());
  498.     }
  500.     // print out overall self-hit rate
  501.     int nSelfHits = 0;
  502.     static const double threshold = 0.01;
  503.     for (row=0; row<multiple->NRows(); ++row) {
  504.         if (multiple->GetRowDouble(row) >= 0.0 && multiple->GetRowDouble(row) <= threshold)
  505.             ++nSelfHits;
  506.     }
  507.     INFOMSG("Self hits with E-value <= " << setprecision(3) << threshold << ": "
  508.         << (100.0*nSelfHits/multiple->NRows()) << "% ("
  509.         << nSelfHits << '/' << multiple->NRows() << ')' << setprecision(6));
  511.     BLASTOptionDelete(options);
  512.     masterSeqInt->id = NULL;    // don't free Seq-id, since it belongs to the Bioseq
  513.     SeqIntFree(masterSeqInt);
  514.     MemFree(masterSeqLoc);
  515.     slaveSeqInt->id = NULL;
  516.     SeqIntFree(slaveSeqInt);
  517.     MemFree(slaveSeqLoc);
  518. }
  520. END_SCOPE(Cn3D)
  522. /*
  523. * ---------------------------------------------------------------------------
  524. * $Log: cn3d_blast.cpp,v $
  525. * Revision 1000.2  2004/06/01 18:28:09  gouriano
  526. * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.34
  527. *
  528. * Revision 1.34  2004/05/21 21:41:39  gorelenk
  529. * Added PCH ncbi_pch.hpp
  530. *
  531. * Revision 1.33  2004/03/15 18:19:23  thiessen
  532. * prefer prefix vs. postfix ++/-- operators
  533. *
  534. * Revision 1.32  2004/02/19 17:04:49  thiessen
  535. * remove cn3d/ from include paths; add pragma to disable annoying msvc warning
  536. *
  537. * Revision 1.31  2003/07/14 18:37:07  thiessen
  538. * change GetUngappedAlignedBlocks() param types; other syntax changes
  539. *
  540. * Revision 1.30  2003/05/29 16:38:27  thiessen
  541. * set db length for CDD 1.62
  542. *
  543. * Revision 1.29  2003/02/03 19:20:02  thiessen
  544. * format changes: move CVS Log to bottom of file, remove std:: from .cpp files, and use new diagnostic macros
  545. *
  546. * Revision 1.28  2003/01/31 17:18:58  thiessen
  547. * many small additions and changes...
  548. *
  549. * Revision 1.27  2003/01/23 20:03:05  thiessen
  550. * add BLAST Neighbor algorithm
  551. *
  552. * Revision 1.26  2003/01/22 14:47:30  thiessen
  553. * cache PSSM in BlockMultipleAlignment
  554. *
  555. * Revision 1.25  2002/12/20 02:51:46  thiessen
  556. * fix Prinf to self problems
  557. *
  558. * Revision 1.24  2002/10/25 19:00:02  thiessen
  559. * retrieve VAST alignment from vastalign.cgi on structure import
  560. *
  561. * Revision 1.23  2002/09/19 14:21:03  thiessen
  562. * add search space size calculation to match scores with rpsblast (finally!)
  563. *
  564. * Revision 1.22  2002/09/05 17:39:53  thiessen
  565. * add bit score printout for BLAST/PSSM
  566. *
  567. * Revision 1.21  2002/09/05 13:04:33  thiessen
  568. * restore output stream precision
  569. *
  570. * Revision 1.20  2002/09/04 15:51:52  thiessen
  571. * turn off options->tweak_parameters
  572. *
  573. * Revision 1.19  2002/08/30 16:52:10  thiessen
  574. * progress on trying to match scores with RPS-BLAST
  575. *
  576. * Revision 1.18  2002/08/15 22:13:13  thiessen
  577. * update for wx2.3.2+ only; add structure pick dialog; fix MultitextDialog bug
  578. *
  579. * Revision 1.17  2002/08/01 12:51:36  thiessen
  580. * add E-value display to block aligner
  581. *
  582. * Revision 1.16  2002/07/26 15:28:45  thiessen
  583. * add Alejandro's block alignment algorithm
  584. *
  585. * Revision 1.15  2002/07/23 15:46:49  thiessen
  586. * print out more BLAST info; tweak save file name
  587. *
  588. * Revision 1.14  2002/07/12 13:24:10  thiessen
  589. * fixes for PSSM creation to agree with cddumper/RPSBLAST
  590. *
  591. * Revision 1.13  2002/06/17 19:31:54  thiessen
  592. * set db_length in blast options
  593. *
  594. * Revision 1.12  2002/06/13 14:54:07  thiessen
  595. * add sort by self-hit
  596. *
  597. * Revision 1.11  2002/06/13 13:32:39  thiessen
  598. * add self-hit calculation
  599. *
  600. * Revision 1.10  2002/05/22 17:17:08  thiessen
  601. * progress on BLAST interface ; change custom spin ctrl implementation
  602. *
  603. * Revision 1.9  2002/05/17 19:10:27  thiessen
  604. * preliminary range restriction for BLAST/PSSM
  605. *
  606. * Revision 1.8  2002/05/07 20:22:47  thiessen
  607. * fix for BLAST/PSSM
  608. *
  609. * Revision 1.7  2002/05/02 18:40:25  thiessen
  610. * do BLAST/PSSM for debug builds only, for testing
  611. *
  612. * Revision 1.6  2002/03/28 14:06:02  thiessen
  613. * preliminary BLAST/PSSM ; new CD startup style
  614. *
  615. * Revision 1.5  2001/11/27 16:26:07  thiessen
  616. * major update to data management system
  617. *
  618. * Revision 1.4  2001/10/18 19:57:32  thiessen
  619. * fix redundant creation of C bioseqs
  620. *
  621. * Revision 1.3  2001/09/27 15:37:58  thiessen
  622. * decouple sequence import and BLAST
  623. *
  624. * Revision 1.2  2001/09/19 22:55:39  thiessen
  625. * add preliminary net import and BLAST
  626. *
  627. * Revision 1.1  2001/09/18 03:10:45  thiessen
  628. * add preliminary sequence import pipeline
  629. *
  630. */