cn3d_threader.cpp

ncbi源码

    if (gibScd) FreeGibScd(gibScd);
    if (trajectory) delete[] trajectory;

    return retval;
}

static double CalculatePSSMScore(const BlockMultipleAlignment::UngappedAlignedBlockList& aBlocks,
    int row, const vector < int >& residueNumbers, const Seq_Mtf *seqMtf)
{
    double score = 0.0;
    BlockMultipleAlignment::UngappedAlignedBlockList::const_iterator b, be = aBlocks.end();
    const Block::Range *masterRange, *slaveRange;
    int i;

    for (b=aBlocks.begin(); b!=be; ++b) {
        masterRange = (*b)->GetRangeOfRow(0);
        slaveRange = (*b)->GetRangeOfRow(row);
        for (i=0; i<(*b)->width; ++i)
            if (residueNumbers[slaveRange->from + i] >= 0)
                score += seqMtf->ww[masterRange->from + i][residueNumbers[slaveRange->from + i]];
    }

//    TESTMSG("PSSM score for row " << row << ": " << score);
    return score;
}

static double CalculateContactScore(const BlockMultipleAlignment *multiple,
    int row, const vector < int >& residueNumbers, const Fld_Mtf *fldMtf, const Rcx_Ptl *rcxPtl)
{
    double score = 0.0;
    int i, seqIndex1, seqIndex2, resNum1, resNum2, dist;

    // for each res-res contact, convert seqIndexes of master into corresponding seqIndexes
    // of slave if they're aligned; add contact energies if so
    for (i=0; i<fldMtf->rrc.n; ++i) {
        seqIndex1 = multiple->GetAlignedSlaveIndex(fldMtf->rrc.r1[i], row);
        if (seqIndex1 < 0) continue;
        seqIndex2 = multiple->GetAlignedSlaveIndex(fldMtf->rrc.r2[i], row);
        if (seqIndex2 < 0) continue;

        resNum1 = residueNumbers[seqIndex1];
        resNum2 = residueNumbers[seqIndex2];
        if (resNum1 < 0 || resNum2 < 0) continue;

        dist = fldMtf->rrc.d[i];
        score += rcxPtl->rre[dist][resNum1][resNum2] + rcxPtl->re[dist][resNum1] + rcxPtl->re[dist][resNum2];
    }

    // ditto for res-pep contacts - except only one slave residue to look up; 2nd is always peptide group
    for (i=0; i<fldMtf->rpc.n; ++i) {
        seqIndex1 = multiple->GetAlignedSlaveIndex(fldMtf->rpc.r1[i], row);
        if (seqIndex1 < 0) continue;

        // peptides are only counted if both contributing master residues are aligned
        if (fldMtf->rpc.p2[i] >= multiple->GetMaster()->Length() - 1 ||
            !multiple->IsAligned(0, fldMtf->rpc.p2[i]) ||
            !multiple->IsAligned(0, fldMtf->rpc.p2[i] + 1)) continue;

        resNum1 = residueNumbers[seqIndex1];
        if (resNum1 < 0) continue;
        resNum2 = NUM_RES_TYPES - 1; // peptide group

        dist = fldMtf->rpc.d[i];
        score += rcxPtl->rre[dist][resNum1][resNum2] + rcxPtl->re[dist][resNum1] + rcxPtl->re[dist][resNum2];
    }

//    TESTMSG("Contact score for row " << row << ": " << score);
    return score;
}

bool Threader::CalculateScores(const BlockMultipleAlignment *multiple, double weightPSSM)
{
    Seq_Mtf *seqMtf = NULL;
    Rcx_Ptl *rcxPtl = NULL;
    Fld_Mtf *fldMtf = NULL;
    BlockMultipleAlignment::UngappedAlignedBlockList aBlocks;
    vector < int > residueNumbers;
    bool retval = false;
    int row;

    // create contact lists
    if (weightPSSM < 1.0 && (!multiple->GetMaster()->molecule ||
            multiple->GetMaster()->molecule->parentSet->isAlphaOnly)) {
        ERRORMSG("Can't use contact potential on non-structured master, or alpha-only (virtual bond) models!");
        goto cleanup;
    }
    if (weightPSSM < 1.0 && !(fldMtf = CreateFldMtf(multiple->GetMaster()))) goto cleanup;

    // create PSSM
    if (weightPSSM > 0.0 && !(seqMtf = CreateSeqMtf(multiple, weightPSSM, NULL))) goto cleanup;

    // create potential
    if (weightPSSM < 1.0 && !(rcxPtl = CreateRcxPtl(1.0 - weightPSSM))) goto cleanup;

    // get aligned blocks
    multiple->GetUngappedAlignedBlocks(&aBlocks);

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

        // get sequence's residue numbers
        const Sequence *seq = multiple->GetSequenceOfRow(row);
        residueNumbers.resize(seq->Length());
        for (int i=0; i<seq->Length(); ++i)
            residueNumbers[i] = LookupThreaderResidueNumberFromCharacterAbbrev(seq->sequenceString[i]);

        // sum score types (weightPSSM already built into seqMtf & rcxPtl)
        double
            scorePSSM = (weightPSSM > 0.0) ?
                CalculatePSSMScore(aBlocks, row, residueNumbers, seqMtf) : 0.0,
            scoreContacts = (weightPSSM < 1.0) ?
                CalculateContactScore(multiple, row, residueNumbers, fldMtf, rcxPtl) : 0.0,
            score = (scorePSSM + scoreContacts) / SCALING_FACTOR;

        // set score in alignment rows (for sorting and status line display)
        multiple->SetRowDouble(row, score);
        CNcbiOstrstream oss;
        oss << "PSSM+Contact score (PSSM x" << weightPSSM << "): " << score << '\0';
        multiple->SetRowStatusLine(row, oss.str());
        delete oss.str();
    }

    retval = true;

cleanup:
    if (seqMtf) FreeSeqMtf(seqMtf);
    if (rcxPtl) FreeRcxPtl(rcxPtl);
    return retval;
}

int Threader::GetGeometryViolations(const BlockMultipleAlignment *multiple,
    GeometryViolationsForRow *violations)
{
    Fld_Mtf *fldMtf = NULL;

    // create contact lists
    if (!multiple->GetMaster()->molecule || multiple->GetMaster()->molecule->parentSet->isAlphaOnly) {
        ERRORMSG("Can't use contact potential on non-structured master, or alpha-only (virtual bond) models!");
        return false;
    }
    if (!(fldMtf = CreateFldMtf(multiple->GetMaster()))) return false;

    violations->clear();
    violations->resize(multiple->NRows());

    // look for too-short regions between aligned blocks
    BlockMultipleAlignment::UngappedAlignedBlockList aBlocks;
    multiple->GetUngappedAlignedBlocks(&aBlocks);
    BlockMultipleAlignment::UngappedAlignedBlockList::const_iterator b, be = aBlocks.end(), n;
    int nViolations = 0;
    const Block::Range *thisRange, *nextRange, *thisMaster, *nextMaster;
    for (b=aBlocks.begin(); b!=be; ++b) {
        n = b;
        ++n;
        if (n == be) break;
        thisMaster = (*b)->GetRangeOfRow(0);
        nextMaster = (*n)->GetRangeOfRow(0);

        for (int row=1; row<multiple->NRows(); ++row) {
            thisRange = (*b)->GetRangeOfRow(row);
            nextRange = (*n)->GetRangeOfRow(row);

            // violation found
            if (nextRange->from - thisRange->to - 1 < fldMtf->mll[nextMaster->from][thisMaster->to]) {
                (*violations)[row].push_back(make_pair(thisRange->to, nextRange->from));
                ++nViolations;
            }
        }
    }

//    TESTMSG("Found " << nViolations << " geometry violations");
    return nViolations;
}

int Threader::EstimateNRandomStarts(const BlockMultipleAlignment *coreAlignment,
    const BlockMultipleAlignment *toBeThreaded)
{
    int nBlocksToAlign = 0;
    BlockMultipleAlignment::UngappedAlignedBlockList multipleABlocks, pairwiseABlocks;
    coreAlignment->GetUngappedAlignedBlocks(&multipleABlocks);
    toBeThreaded->GetUngappedAlignedBlocks(&pairwiseABlocks);

    // if a block in the multiple is *not* contained in the pairwise (looking at master coords),
    // then it'll probably be realigned upon threading
    BlockMultipleAlignment::UngappedAlignedBlockList::const_iterator
        m, me = multipleABlocks.end(), p, pe = pairwiseABlocks.end();
    const Block::Range *multipleRange, *pairwiseRange;
    for (m=multipleABlocks.begin(); m!=me; ++m) {
        multipleRange = (*m)->GetRangeOfRow(0);
        bool realignBlock = true;
        for (p=pairwiseABlocks.begin(); p!=pe; ++p) {
            pairwiseRange = (*p)->GetRangeOfRow(0);
            if (pairwiseRange->from <= multipleRange->from && pairwiseRange->to >= multipleRange->to) {
                realignBlock = false;
                break;
            }
        }
        if (realignBlock) ++nBlocksToAlign;
    }

    if (nBlocksToAlign <= 1)
        return 1;
    else
        // round to nearest integer
        return (int) (exp(1.5 + 0.25432 * nBlocksToAlign) + 0.5);
}

END_SCOPE(Cn3D)

/*
* ---------------------------------------------------------------------------
* $Log: cn3d_threader.cpp,v $
* Revision 1000.3  2004/06/01 18:28:21  gouriano
* PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.43
*
* Revision 1.43  2004/05/21 21:41:39  gorelenk
* Added PCH ncbi_pch.hpp
*
* Revision 1.42  2004/03/15 17:59:20  thiessen
* prefer prefix vs. postfix ++/-- operators
*
* Revision 1.41  2004/03/15 17:33:12  thiessen
* prefer prefix vs. postfix ++/-- operators
*
* Revision 1.40  2004/02/19 17:04:51  thiessen
* remove cn3d/ from include paths; add pragma to disable annoying msvc warning
*
* Revision 1.39  2003/11/04 18:09:17  thiessen
* rearrange headers for OSX build
*
* Revision 1.38  2003/07/14 18:37:07  thiessen
* change GetUngappedAlignedBlocks() param types; other syntax changes
*
* Revision 1.37  2003/06/18 11:38:42  thiessen
* add another trace message
*
* Revision 1.36  2003/04/14 18:13:58  thiessen
* retry pseudocounts until matrix is constructed okay
*
* Revision 1.35  2003/02/03 19:20:03  thiessen
* format changes: move CVS Log to bottom of file, remove std:: from .cpp files, and use new diagnostic macros
*
* Revision 1.34  2003/01/23 20:03:05  thiessen
* add BLAST Neighbor algorithm
*
* Revision 1.33  2002/10/08 12:35:42  thiessen
* use delete[] for arrays
*
* Revision 1.32  2002/08/15 22:13:14  thiessen
* update for wx2.3.2+ only; add structure pick dialog; fix MultitextDialog bug
*
* Revision 1.31  2002/07/26 15:28:47  thiessen
* add Alejandro's block alignment algorithm
*
* Revision 1.30  2002/07/12 13:24:10  thiessen
* fixes for PSSM creation to agree with cddumper/RPSBLAST
*
* Revision 1.29  2002/05/26 21:58:46  thiessen
* add CddDegapSeqAlign to PSSM generator
*
* Revision 1.28  2002/03/28 14:06:02  thiessen
* preliminary BLAST/PSSM ; new CD startup style
*
* Revision 1.27  2002/03/19 18:47:57  thiessen
* small bug fixes; remember PSSM weight
*
* Revision 1.26  2002/02/21 22:01:49  thiessen
* remember alignment range on demotion
*
* Revision 1.25  2002/02/21 12:26:29  thiessen
* fix row delete bug ; remember threader options
*
* Revision 1.24  2001/10/08 00:00:09  thiessen
* estimate threader N random starts; edit CDD name
*
* Revision 1.23  2001/09/27 15:37:58  thiessen
* decouple sequence import and BLAST
*
* Revision 1.22  2001/06/21 02:02:33  thiessen
* major update to molecule identification and highlighting ; add toggle highlight (via alt)
*
* Revision 1.21  2001/06/02 17:22:45  thiessen
* fixes for GCC
*
* Revision 1.20  2001/06/01 21:48:26  thiessen
* add terminal cutoff to threading
*
* Revision 1.19  2001/05/31 18:47:07  thiessen
* add preliminary style dialog; remove LIST_TYPE; add thread single and delete all; misc tweaks
*
* Revision 1.18  2001/05/24 21:38:41  thiessen
* fix threader options initial values
*
* Revision 1.17  2001/05/15 23:48:36  thiessen
* minor adjustments to compile under Solaris/wxGTK
*
* Revision 1.16  2001/05/15 14:57:55  thiessen
* add cn3d_tools; bring up log window when threading starts
*
* Revision 1.15  2001/05/11 13:45:06  thiessen
* set up data directory
*
* Revision 1.14  2001/05/11 02:10:42  thiessen
* add better merge fail indicators; tweaks to windowing/taskbar
*
* Revision 1.13  2001/04/13 18:50:54  thiessen
* fix for when threader returns fewer than requested results
*
* Revision 1.12  2001/04/12 18:54:39  thiessen
* fix memory leak for PSSM-only threading
*
* Revision 1.11  2001/04/12 18:10:00  thiessen
* add block freezing
*
* Revision 1.10  2001/04/05 22:55:35  thiessen
* change bg color handling ; show geometry violations
*
* Revision 1.9  2001/04/04 00:27:14  thiessen
* major update - add merging, threader GUI controls
*
* Revision 1.8  2001/03/30 14:43:41  thiessen
* show threader scores in status line; misc UI tweaks
*
* Revision 1.7  2001/03/30 03:07:34  thiessen
* add threader score calculation & sorting
*
* Revision 1.6  2001/03/29 15:35:55  thiessen
* remove GetAtom warnings
*
* Revision 1.5  2001/03/28 23:02:16  thiessen
* first working full threading
*
* Revision 1.4  2001/03/23 23:31:56  thiessen
* keep atom info around even if coords not all present; mainly for disulfide parsing in virtual models
*
* Revision 1.3  2001/03/22 00:33:16  thiessen
* initial threading working (PSSM only); free color storage in undo stack
*
* Revision 1.2  2001/02/13 01:03:56  thiessen
* backward-compatible domain ID's in output; add ability to delete rows
*
* Revision 1.1  2001/02/08 23:01:49  thiessen
* hook up C-toolkit stuff for threading; working PSSM calculation
*
*/