cn3d_blast.cpp

ncbi源码

            TRACEMSG("Calculated raw score by BLOSUM62: " << raw_bl62);
#endif

            // get scores
            if (sa.IsSetScore()) {
                wxString scores;
                scores.Printf("BLAST result for %s vs. %s:",
                    masterSeq->identifier->ToString().c_str(),
                    slaveSeq->identifier->ToString().c_str());

                CSeq_align::TScore::const_iterator sc, sce = sa.GetScore().end();
                for (sc=sa.GetScore().begin(); sc!=sce; ++sc) {
                    if ((*sc)->IsSetId() && (*sc)->GetId().IsStr()) {

                        // E-value (put in status line and double values)
                        if ((*sc)->GetValue().IsReal() && (*sc)->GetId().GetStr() == "e_value") {
                            newAlignment->SetRowDouble(0, (*sc)->GetValue().GetReal());
                            newAlignment->SetRowDouble(1, (*sc)->GetValue().GetReal());
                            wxString status;
                            status.Printf("E-value %g", (*sc)->GetValue().GetReal());
                            newAlignment->SetRowStatusLine(0, status.c_str());
                            newAlignment->SetRowStatusLine(1, status.c_str());
                            wxString tmp = scores;
                            scores.Printf("%s E-value: %g", tmp.c_str(), (*sc)->GetValue().GetReal());
                        }

                        // raw score
                        if ((*sc)->GetValue().IsInt() && (*sc)->GetId().GetStr() == "score") {
                            wxString tmp = scores;
                            scores.Printf("%s raw: %i", tmp.c_str(), (*sc)->GetValue().GetInt());
                        }

                        // bit score
                        if ((*sc)->GetValue().IsReal() && (*sc)->GetId().GetStr() == "bit_score") {
                            wxString tmp = scores;
                            scores.Printf("%s bit score: %g", tmp.c_str(), (*sc)->GetValue().GetReal());
                        }
                    }
                }

                INFOMSG(scores.c_str());
            }
        }

        // finalize and and add new alignment to list
        if (newAlignment->AddUnalignedBlocks() && newAlignment->UpdateBlockMapAndColors(false)) {
            newAlignments->push_back(newAlignment);
        } else {
            ERRORMSG("error finalizing alignment");
            delete newAlignment;
        }
    }

    BLASTOptionDelete(options);
    masterSeqInt->id = NULL;    // don't free Seq-id, since it belongs to the Bioseq
    SeqIntFree(masterSeqInt);
    MemFree(masterSeqLoc);
    slaveSeqInt->id = NULL;
    SeqIntFree(slaveSeqInt);
    MemFree(slaveSeqLoc);
}

void BLASTer::CalculateSelfHitScores(const BlockMultipleAlignment *multiple)
{
    if (!multiple) {
        ERRORMSG("NULL multiple alignment");
        return;
    }

    // set up BLAST stuff
    BLAST_OptionsBlkPtr options = CreateBlastOptionsBlk();
    if (!options) {
        ERRORMSG("BLASTOptionNew() failed");
        return;
    }

    // get master sequence
    const Sequence *masterSeq = multiple->GetMaster();
    Bioseq *masterBioseq = masterSeq->parentSet->GetOrCreateBioseq(masterSeq);

    // create Seq-loc interval for master and slave
    SeqLocPtr masterSeqLoc = (SeqLocPtr) MemNew(sizeof(SeqLoc));
    masterSeqLoc->choice = SEQLOC_INT;
    SeqIntPtr masterSeqInt = SeqIntNew();
    masterSeqLoc->data.ptrvalue = masterSeqInt;
    masterSeqInt->id = masterBioseq->id;
    BlockMultipleAlignment::UngappedAlignedBlockList uaBlocks;
    multiple->GetUngappedAlignedBlocks(&uaBlocks);
    if (uaBlocks.size() == 0) {
        ERRORMSG("Self-hit requires at least one aligned block");
        return;
    }
    int excess = 0;
    if (!RegistryGetInteger(REG_ADVANCED_SECTION, REG_FOOTPRINT_RES, &excess))
        WARNINGMSG("Can't get footprint excess residues from registry");
    masterSeqInt->from = uaBlocks.front()->GetRangeOfRow(0)->from - excess;
    if (masterSeqInt->from < 0) masterSeqInt->from = 0;
    masterSeqInt->to = uaBlocks.back()->GetRangeOfRow(0)->to + excess;
    if (masterSeqInt->to >= masterSeq->Length()) masterSeqInt->to = masterSeq->Length() - 1;

    SeqLocPtr slaveSeqLoc = (SeqLocPtr) MemNew(sizeof(SeqLoc));
    slaveSeqLoc->choice = SEQLOC_INT;
    SeqIntPtr slaveSeqInt = SeqIntNew();
    slaveSeqLoc->data.ptrvalue = slaveSeqInt;

    // create BLAST_Matrix if using PSSM from multiple
    const BLAST_Matrix *BLASTmatrix = multiple->GetPSSM();

    string err;
    int row;
    for (row=0; row<multiple->NRows(); ++row) {

        // get C Bioseq for each slave
        const Sequence *slaveSeq = multiple->GetSequenceOfRow(row);
        Bioseq *slaveBioseq = slaveSeq->parentSet->GetOrCreateBioseq(slaveSeq);

        // setup Seq-loc interval for slave
        slaveSeqInt->id = slaveBioseq->id;
        slaveSeqInt->from = uaBlocks.front()->GetRangeOfRow(row)->from - excess;
        if (slaveSeqInt->from < 0) slaveSeqInt->from = 0;
        slaveSeqInt->to = uaBlocks.back()->GetRangeOfRow(row)->to + excess;
        if (slaveSeqInt->to >= slaveSeq->Length()) slaveSeqInt->to = slaveSeq->Length() - 1;
//        TESTMSG("for BLAST: master range " <<
//                (masterSeqInt->from + 1) << " to " << (masterSeqInt->to + 1) << ", slave range " <<
//                (slaveSeqInt->from + 1) << " to " << (slaveSeqInt->to + 1));

        // set search space size
        SetEffectiveSearchSpaceSize(options, slaveSeqInt->to - slaveSeqInt->from + 1);

        // actually do the BLAST alignment

//        TRACEMSG("calling BlastTwoSequencesByLocWithCallback()");
        SeqAlign *salp = BlastTwoSequencesByLocWithCallback(
            masterSeqLoc, slaveSeqLoc,
            "blastp", options,
            NULL, NULL, NULL,
            const_cast<BLAST_Matrix*>(BLASTmatrix));

//        TRACEMSG("calling B2SPssmMultipleQueries()");
//        SeqLocPtr target[1];
//        target[0] = slaveSeqLoc;
//        options->searchsp_eff = 0;
//        SeqAlignPtr *psalp = B2SPssmMultipleQueries(masterSeqLoc,
//            const_cast<BLAST_Matrix*>(BLASTmatrix), target, 1, options);
//        SeqAlign *salp = psalp[0];

        // process the result
        double score = -1.0;
        if (salp) {
            // convert C SeqAlign to C++ for convenience
            CSeq_align sa;
            bool okay = ConvertAsnFromCToCPP(salp, (AsnWriteFunc) SeqAlignAsnWrite, &sa, &err);
            SeqAlignSetFree(salp);
            if (!okay) {
                ERRORMSG("Conversion of SeqAlign to C++ object failed");
                continue;
            }

            // get score
            if (sa.IsSetScore()) {
                CSeq_align::TScore::const_iterator sc, sce = sa.GetScore().end();
                for (sc=sa.GetScore().begin(); sc!=sce; ++sc) {
                    if ((*sc)->GetValue().IsReal() && (*sc)->IsSetId() &&
                        (*sc)->GetId().IsStr() && (*sc)->GetId().GetStr() == "e_value") {
                        score = (*sc)->GetValue().GetReal();
                        break;
                    }
                }
            }
            if (score < 0.0) ERRORMSG("Got back Seq-align with no E-value");

        }

        // set score in row
        multiple->SetRowDouble(row, score);
        wxString status;
        if (score >= 0.0)
            status.Printf("Self hit E-value: %g", score);
        else
            status = "No detectable self hit";
        multiple->SetRowStatusLine(row, status.c_str());
    }

    // print out overall self-hit rate
    int nSelfHits = 0;
    static const double threshold = 0.01;
    for (row=0; row<multiple->NRows(); ++row) {
        if (multiple->GetRowDouble(row) >= 0.0 && multiple->GetRowDouble(row) <= threshold)
            ++nSelfHits;
    }
    INFOMSG("Self hits with E-value <= " << setprecision(3) << threshold << ": "
        << (100.0*nSelfHits/multiple->NRows()) << "% ("
        << nSelfHits << '/' << multiple->NRows() << ')' << setprecision(6));

    BLASTOptionDelete(options);
    masterSeqInt->id = NULL;    // don't free Seq-id, since it belongs to the Bioseq
    SeqIntFree(masterSeqInt);
    MemFree(masterSeqLoc);
    slaveSeqInt->id = NULL;
    SeqIntFree(slaveSeqInt);
    MemFree(slaveSeqLoc);
}

END_SCOPE(Cn3D)

/*
* ---------------------------------------------------------------------------
* $Log: cn3d_blast.cpp,v $
* Revision 1000.2  2004/06/01 18:28:09  gouriano
* PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.34
*
* Revision 1.34  2004/05/21 21:41:39  gorelenk
* Added PCH ncbi_pch.hpp
*
* Revision 1.33  2004/03/15 18:19:23  thiessen
* prefer prefix vs. postfix ++/-- operators
*
* Revision 1.32  2004/02/19 17:04:49  thiessen
* remove cn3d/ from include paths; add pragma to disable annoying msvc warning
*
* Revision 1.31  2003/07/14 18:37:07  thiessen
* change GetUngappedAlignedBlocks() param types; other syntax changes
*
* Revision 1.30  2003/05/29 16:38:27  thiessen
* set db length for CDD 1.62
*
* Revision 1.29  2003/02/03 19:20:02  thiessen
* format changes: move CVS Log to bottom of file, remove std:: from .cpp files, and use new diagnostic macros
*
* Revision 1.28  2003/01/31 17:18:58  thiessen
* many small additions and changes...
*
* Revision 1.27  2003/01/23 20:03:05  thiessen
* add BLAST Neighbor algorithm
*
* Revision 1.26  2003/01/22 14:47:30  thiessen
* cache PSSM in BlockMultipleAlignment
*
* Revision 1.25  2002/12/20 02:51:46  thiessen
* fix Prinf to self problems
*
* Revision 1.24  2002/10/25 19:00:02  thiessen
* retrieve VAST alignment from vastalign.cgi on structure import
*
* Revision 1.23  2002/09/19 14:21:03  thiessen
* add search space size calculation to match scores with rpsblast (finally!)
*
* Revision 1.22  2002/09/05 17:39:53  thiessen
* add bit score printout for BLAST/PSSM
*
* Revision 1.21  2002/09/05 13:04:33  thiessen
* restore output stream precision
*
* Revision 1.20  2002/09/04 15:51:52  thiessen
* turn off options->tweak_parameters
*
* Revision 1.19  2002/08/30 16:52:10  thiessen
* progress on trying to match scores with RPS-BLAST
*
* Revision 1.18  2002/08/15 22:13:13  thiessen
* update for wx2.3.2+ only; add structure pick dialog; fix MultitextDialog bug
*
* Revision 1.17  2002/08/01 12:51:36  thiessen
* add E-value display to block aligner
*
* Revision 1.16  2002/07/26 15:28:45  thiessen
* add Alejandro's block alignment algorithm
*
* Revision 1.15  2002/07/23 15:46:49  thiessen
* print out more BLAST info; tweak save file name
*
* Revision 1.14  2002/07/12 13:24:10  thiessen
* fixes for PSSM creation to agree with cddumper/RPSBLAST
*
* Revision 1.13  2002/06/17 19:31:54  thiessen
* set db_length in blast options
*
* Revision 1.12  2002/06/13 14:54:07  thiessen
* add sort by self-hit
*
* Revision 1.11  2002/06/13 13:32:39  thiessen
* add self-hit calculation
*
* Revision 1.10  2002/05/22 17:17:08  thiessen
* progress on BLAST interface ; change custom spin ctrl implementation
*
* Revision 1.9  2002/05/17 19:10:27  thiessen
* preliminary range restriction for BLAST/PSSM
*
* Revision 1.8  2002/05/07 20:22:47  thiessen
* fix for BLAST/PSSM
*
* Revision 1.7  2002/05/02 18:40:25  thiessen
* do BLAST/PSSM for debug builds only, for testing
*
* Revision 1.6  2002/03/28 14:06:02  thiessen
* preliminary BLAST/PSSM ; new CD startup style
*
* Revision 1.5  2001/11/27 16:26:07  thiessen
* major update to data management system
*
* Revision 1.4  2001/10/18 19:57:32  thiessen
* fix redundant creation of C bioseqs
*
* Revision 1.3  2001/09/27 15:37:58  thiessen
* decouple sequence import and BLAST
*
* Revision 1.2  2001/09/19 22:55:39  thiessen
* add preliminary net import and BLAST
*
* Revision 1.1  2001/09/18 03:10:45  thiessen
* add preliminary sequence import pipeline
*
*/