cn3d_threader.cpp

ncbi源码

    }
}

static void TranslateContacts(const Threader::ContactList& resResContacts,
    const Threader::ContactList& resPepContacts, Fld_Mtf *fldMtf)
{
    int i;
    Threader::ContactList::const_iterator c;
    for (i=0, c=resResContacts.begin(); i<resResContacts.size(); ++i, ++c) {
        fldMtf->rrc.r1[i] = c->vc1 / 2;  // threader coord points to (res,pep) pair
        fldMtf->rrc.r2[i] = c->vc2 / 2;
        fldMtf->rrc.d[i] = c->distanceBin;
    }
    for (i=0, c=resPepContacts.begin(); i<resPepContacts.size(); ++i, ++c) {
        fldMtf->rpc.r1[i] = c->vc1 / 2;
        fldMtf->rpc.p2[i] = c->vc2 / 2;
        fldMtf->rpc.d[i] = c->distanceBin;
    }
}

// for sorting contacts
inline bool operator < (const Threader::Contact& c1, const Threader::Contact& c2)
{
    return (c1.vc1 < c2.vc1 || (c1.vc1 == c2.vc1 && c1.vc2 < c2.vc2));
}

static void GetMinimumLoopLengths(const Molecule *mol, const AtomSet *atomSet, Fld_Mtf *fldMtf)
{
    int i, j;
    const AtomCoord *a1, *a2;
    Molecule::ResidueMap::const_iterator r1, r2, re = mol->residues.end();
    for (r1=mol->residues.begin(), i=0; r1!=re; ++r1, ++i) {

        if (r1->second->alphaID == Residue::NO_ALPHA_ID)
            a1 = NULL;
        else {
            AtomPntr ap1(mol->id, r1->second->id, r1->second->alphaID);
            a1 = atomSet->GetAtom(ap1, true, true);   // 'true' means just use first alt coord
        }

        for (r2=r1, j=i; r2!=re; ++r2, ++j) {

            if (i == j) {
                fldMtf->mll[i][j] = 0;
            } else {
                if (r2->second->alphaID == Residue::NO_ALPHA_ID)
                    a2 = NULL;
                else {
                    AtomPntr ap2(mol->id, r2->second->id, r2->second->alphaID);
                    a2 = atomSet->GetAtom(ap2, true, true);
                }
                fldMtf->mll[i][j] = fldMtf->mll[j][i] =
                    (!a1 || !a2) ? 0 : (int) (((a2->site - a1->site).length() - 2.7) / 3.4);
            }
        }
    }
}

Fld_Mtf * Threader::CreateFldMtf(const Sequence *masterSequence)
{
    if (!masterSequence) return NULL;

    const Molecule *mol = masterSequence->molecule;

    // return cached copy if we've already constructed a Fld_Mtf for this master
    ContactMap::iterator c = mol ? contacts.find(mol) : contacts.find(masterSequence);
    if (c != contacts.end()) return c->second;

    // work-around to allow PSSM-only threading when master has no structure (or only C-alphas)
    Fld_Mtf *fldMtf;
    if (!mol || mol->parentSet->isAlphaOnly) {
        fldMtf = NewFldMtf(masterSequence->Length(), 0, 0);
        contacts[masterSequence] = fldMtf;
        return fldMtf;
    }

    // for convenience so subroutines don't have to keep looking this up... Use first
    // CoordSet if multiple model (e.g., NMR)
    const StructureObject *object;
    if (!mol->GetParentOfType(&object)) return NULL;
    const AtomSet *atomSet = object->coordSets.front()->atomSet;

    // get virtual coordinates for this chain
    VirtualCoordinateList virtualCoordinates;
    GetVirtualCoordinates(mol, atomSet, &virtualCoordinates);

    // check for contacts of virtual coords separated by >= 10 virtual bonds
    ContactList resResContacts, resPepContacts;
    GetContacts(virtualCoordinates, &resResContacts, &resPepContacts);

    // create Fld_Mtf, and store contacts in it
    fldMtf = NewFldMtf(mol->residues.size(), resResContacts.size(), resPepContacts.size());
    resPepContacts.sort();  // not really necessary, but makes same order as Cn3D for comparison/testing
    TranslateContacts(resResContacts, resPepContacts, fldMtf);

    // fill out min. loop lengths
    GetMinimumLoopLengths(mol, atomSet, fldMtf);

    TRACEMSG("created Fld_Mtf for " << mol->identifier->pdbID << " chain '" << (char) mol->identifier->pdbChain << "'");
    contacts[mol] = fldMtf;
    return fldMtf;
}

static BlockMultipleAlignment * CreateAlignmentFromThdTbl(const Thd_Tbl *thdTbl, int nResult,
    const Cor_Def *corDef, BlockMultipleAlignment::SequenceList *sequences, AlignmentManager *alnMgr)
{
    if (corDef->sll.n != thdTbl->nsc || nResult >= thdTbl->n) {
        ERRORMSG("CreateAlignmentFromThdTbl() - inconsistent Thd_Tbl");
        return NULL;
    }

    BlockMultipleAlignment *newAlignment = new BlockMultipleAlignment(sequences, alnMgr);

    // add blocks from threader result
    for (int block=0; block<corDef->sll.n; ++block) {
        UngappedAlignedBlock *aBlock = new UngappedAlignedBlock(newAlignment);
        aBlock->SetRangeOfRow(0,
            corDef->sll.rfpt[block] - thdTbl->no[block][nResult],
            corDef->sll.rfpt[block] + thdTbl->co[block][nResult]);
        aBlock->SetRangeOfRow(1,
            thdTbl->al[block][nResult] - thdTbl->no[block][nResult],
            thdTbl->al[block][nResult] + thdTbl->co[block][nResult]);
        aBlock->width = thdTbl->no[block][nResult] + 1 + thdTbl->co[block][nResult];
        if (!newAlignment->AddAlignedBlockAtEnd(aBlock)) {
            ERRORMSG("CreateAlignmentFromThdTbl() - error adding block");
            delete newAlignment;
            return NULL;
        }
    }

    // finish alignment
    if (!newAlignment->AddUnalignedBlocks() || !newAlignment->UpdateBlockMapAndColors()) {
        ERRORMSG("CreateAlignmentFromThdTbl() - error finishing alignment");
        delete newAlignment;
        return NULL;
    }

    return newAlignment;
}

static bool FreezeIsolatedBlocks(Cor_Def *corDef, const Cor_Def *masterCorDef, const Qry_Seq *qrySeq)
{
    if (!corDef || !masterCorDef || !qrySeq ||
        corDef->sll.n != masterCorDef->sll.n || corDef->sll.n != qrySeq->sac.n) {
        ERRORMSG("FreezeIsolatedBlocks() - bad parameters");
        return false;
    }

    TRACEMSG("freezing blocks...");
    for (int i=0; i<corDef->sll.n; ++i) {

        // default: blocks allowed to grow
        corDef->sll.nomx[i] = masterCorDef->sll.nomx[i];
        corDef->sll.comx[i] = masterCorDef->sll.comx[i];

        // new blocks always allowed to grow
        if (qrySeq->sac.mn[i] < 0 || qrySeq->sac.mx[i] < 0) continue;

        // if an existing block is adjacent to any new (to-be-realigned) block, then allow block's
        // boundaries to grow on that side; otherwise, freeze (isolated) existing block boundaries
        bool adjacentLeft = (i > 0 && (qrySeq->sac.mn[i - 1] < 0 || qrySeq->sac.mx[i - 1] < 0));
        bool adjacentRight = (i < corDef->sll.n - 1 &&
            (qrySeq->sac.mn[i + 1] < 0 || qrySeq->sac.mx[i + 1] < 0));

        if (!adjacentLeft) {
            corDef->sll.nomx[i] = corDef->sll.nomn[i];
//            TESTMSG("block " << i << " fixed N-terminus");
        }
        if (!adjacentRight) {
            corDef->sll.comx[i] = corDef->sll.comn[i];
//            TESTMSG("block " << i << " fixed C-terminus");
        }
    }

    return true;
}

bool Threader::Realign(const ThreaderOptions& options, BlockMultipleAlignment *masterMultiple,
    const AlignmentList *originalAlignments,
    int *nRowsAddedToMultiple, AlignmentList *newAlignments)
{
    *nRowsAddedToMultiple = 0;
    if (!masterMultiple || !originalAlignments || !newAlignments || originalAlignments->size() == 0)
        return false;

    // either calculate no z-scores (0), or calculate z-score for best result (1)
    static const int zscs = 0;

    Seq_Mtf *seqMtf = NULL;
    Cor_Def *corDef = NULL, *masterCorDef = NULL;
    Rcx_Ptl *rcxPtl = NULL;
    Gib_Scd *gibScd = NULL;
    Fld_Mtf *fldMtf = NULL;
    float *trajectory = NULL;
    bool retval = false;

    AlignmentList::const_iterator p, pe = originalAlignments->end();

#ifdef DEBUG_THREADER
    FILE *pFile;
#endif

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

    // create potential and Gibbs schedule
    if (!(rcxPtl = CreateRcxPtl(1.0 - options.weightPSSM))) goto cleanup;
    if (!(gibScd = CreateGibScd(true, options.nRandomStarts))) goto cleanup;
    trajectory = new float[gibScd->ntp];

    // create initial PSSM
    if (!(seqMtf = CreateSeqMtf(masterMultiple, options.weightPSSM, NULL))) goto cleanup;
#ifdef DEBUG_THREADER
    pFile = fopen("Seq_Mtf.debug.txt", "w");
    PrintSeqMtf(seqMtf, pFile);
    fclose(pFile);
#endif

    // create core definition
    if (!(corDef = CreateCorDef(masterMultiple, options.loopLengthMultiplier))) goto cleanup;
    if (options.freezeIsolatedBlocks)   // make a copy to used as an original "master"
        if (!(masterCorDef = CreateCorDef(masterMultiple, options.loopLengthMultiplier))) goto cleanup;

#ifdef DEBUG_THREADER
    pFile = fopen("Fld_Mtf.debug.txt", "w");
    PrintFldMtf(fldMtf, pFile);
    fclose(pFile);
#endif

    for (p=originalAlignments->begin(); p!=pe; ) {

        if ((*p)->NRows() != 2 || (*p)->GetMaster() != masterMultiple->GetMaster()) {
            ERRORMSG("Threader::Realign() - bad pairwise alignment");
            continue;
        }

        Qry_Seq *qrySeq = NULL;
        Thd_Tbl *thdTbl = NULL;
        int success;

        // create query sequence
        if (!(qrySeq = CreateQrySeq(masterMultiple, *p, options.terminalResidueCutoff))) goto cleanup2;
#ifdef DEBUG_THREADER
        pFile = fopen("Qry_Seq.debug.txt", "w");
        PrintQrySeq(qrySeq, pFile);
        fclose(pFile);
#endif

        // freeze block sizes if opted (changes corDef but not masterCorDef or qrySeq)
        if (options.freezeIsolatedBlocks)
            FreezeIsolatedBlocks(corDef, masterCorDef, qrySeq);
#ifdef DEBUG_THREADER
        pFile = fopen("Cor_Def.debug.txt", "w");
        PrintCorDef(corDef, pFile);
        fclose(pFile);
#endif

        // create results storage structure
        thdTbl = NewThdTbl(options.nResultAlignments, corDef->sll.n);

        // actually run the threader (finally!)
        INFOMSG("threading " << (*p)->GetSequenceOfRow(1)->identifier->ToString());
        success = atd(fldMtf, corDef, qrySeq, rcxPtl, gibScd, thdTbl, seqMtf,
            trajectory, zscs, SCALING_FACTOR, (float) options.weightPSSM);

        BlockMultipleAlignment *newAlignment;
        if (success) {
            TRACEMSG("threading succeeded");
#ifdef DEBUG_THREADER
            pFile = fopen("Thd_Tbl.debug.txt", "w");
            PrintThdTbl(thdTbl, pFile);
            fclose(pFile);
#endif
            // create new alignment(s) from threading result; merge or add to list as appropriate
            for (int i=0; i<thdTbl->n; ++i) {

                // skip if this entry is not a real result
                if (thdTbl->tf[i] <= 0) continue;

                BlockMultipleAlignment::SequenceList *sequences = new BlockMultipleAlignment::SequenceList(2);
                sequences->front() = (*p)->GetMaster();
                sequences->back() = (*p)->GetSequenceOfRow(1);
                newAlignment = CreateAlignmentFromThdTbl(thdTbl, i, corDef,
                    sequences, masterMultiple->alignmentManager);
                if (!newAlignment) continue;

                // set scores to show in alignment
                newAlignment->SetRowDouble(0, thdTbl->tg[i]);
                newAlignment->SetRowDouble(1, thdTbl->tg[i]);
                CNcbiOstrstream oss;
                oss << "Threading successful; alignment score before merge: " << thdTbl->tg[i] << '\0';
                newAlignment->SetRowStatusLine(0, oss.str());
                newAlignment->SetRowStatusLine(1, oss.str());
                delete oss.str();

                if (options.mergeAfterEachSequence) {
                    if (masterMultiple->MergeAlignment(newAlignment)) {
                        delete newAlignment; // if merge is successful, we can delete this alignment;
                        ++(*nRowsAddedToMultiple);
                    } else {                 // otherwise keep it
                        newAlignments->push_back(newAlignment);
                    }
                }

                // no merge - add new alignment to list, let calling function deal with it
                else {
                    newAlignments->push_back(newAlignment);
                }
            }
        }

        // threading failed - add old alignment to list so it doesn't get lost
        else {
            TRACEMSG("threading failed!");
            newAlignment = (*p)->Clone();
            newAlignment->SetRowDouble(0, -1.0);
            newAlignment->SetRowDouble(1, -1.0);
            newAlignment->SetRowStatusLine(0, "Threading failed!");
            newAlignment->SetRowStatusLine(1, "Threading failed!");
            newAlignments->push_back(newAlignment);
        }

cleanup2:
        if (qrySeq) FreeQrySeq(qrySeq);
        if (thdTbl) FreeThdTbl(thdTbl);

        ++p;
        if (success && p != pe && options.mergeAfterEachSequence) {
            // re-create PSSM after each merge
            FreeSeqMtf(seqMtf);
            if (!(seqMtf = CreateSeqMtf(masterMultiple, options.weightPSSM, NULL))) goto cleanup;
        }
    }

    retval = true;

cleanup:
    if (seqMtf) FreeSeqMtf(seqMtf);
    if (corDef) FreeCorDef(corDef);
    if (masterCorDef) FreeCorDef(masterCorDef);
    if (rcxPtl) FreeRcxPtl(rcxPtl);