cn3d_threader.cpp

ncbi源码

            ++Count;
            InsideStr = TRUE;
        }
        if (Str[i] == ' ') {
            InsideStr = FALSE;
        }
    }
    return(Count);
}
static void ReadToRowOfEnergies(ifstream& InFile, int NumResTypes) {
    char  Str[1024];
    while (!InFile.eof()) {
        InFile.getline(Str, sizeof(Str));
            if (CountWords(Str) == NumResTypes) {
            break;
        }
    }
}

static const int NUM_RES_TYPES = 21;

Rcx_Ptl * Threader::CreateRcxPtl(double weightContacts)
{
    Rcx_Ptl*  pmf;
    char      *FileName = "ContactPotential";
    char      ResName[32];
    char      Path[512];
    Int4      i, j, k;
    double    temp;

    static const int kNumDistances = 6;
    static const int kPeptideIndex = 20;

    /* open the contact potential for reading */
    StrCpy(Path, GetDataDir().c_str());
    StrCat(Path, FileName);
    auto_ptr<CNcbiIfstream> InFile(new CNcbiIfstream(Path));
    if (!(*InFile)) {
        ERRORMSG("Threader::CreateRcxPtl() - can't open " << Path << " for reading");
        return NULL;
    }

    pmf = NewRcxPtl(NUM_RES_TYPES, kNumDistances, kPeptideIndex);

    /* read in the contact potential */
    for (i=0; i<kNumDistances; ++i) {
        ReadToRowOfEnergies(*InFile, NUM_RES_TYPES);
        if (InFile->eof()) goto error;
        for (j=0; j<NUM_RES_TYPES; ++j) {
            InFile->getline(ResName, sizeof(ResName), ' ');  /* skip residue name */
            if (InFile->eof()) goto error;
            for (k=0; k<NUM_RES_TYPES; ++k) {
                *InFile >> temp;
                if (InFile->eof()) goto error;
                pmf->rre[i][j][k] = ThrdRound(temp*SCALING_FACTOR*weightContacts);
            }
        }
    }

    /* read in the hydrophobic energies */
    ReadToRowOfEnergies(*InFile, kNumDistances);
    for (i=0; i<NUM_RES_TYPES; ++i) {
        InFile->getline(ResName, sizeof(ResName), ' ');  /* skip residue name */
        if (InFile->eof()) goto error;
        for (j=0; j<kNumDistances; ++j) {
            *InFile >> temp;
            if (InFile->eof()) goto error;
            pmf->re[j][i] = ThrdRound(temp*SCALING_FACTOR*weightContacts);
        }
    }

    /* calculate sum of pair energies plus hydrophobic energies */
    for(i=0; i<kNumDistances; ++i) {
        for(j=0; j<NUM_RES_TYPES; ++j) {
            for(k=0; k<NUM_RES_TYPES; ++k) {
                pmf->rrt[i][j][k] = pmf->rre[i][j][k] + pmf->re[i][j] + pmf->re[i][k];
            }
        }
    }

    return(pmf);

error:
    ERRORMSG("Threader::CreateRcxPtl() - error parsing " << FileName);
    FreeRcxPtl(pmf);
    return NULL;
}

/*----------------------------------------------------------------------------
 *  set up the annealing parameters. (more code swiped from DDV)
 *  hard-coded for now.  later we can move these parameters to a file.
 *---------------------------------------------------------------------------*/
Gib_Scd * Threader::CreateGibScd(bool fast, int nRandomStarts)
{
    Gib_Scd*  gsp;
    int       NumTrajectoryPoints;

    static const int kNumTempSteps = 3;

    gsp = NewGibScd(kNumTempSteps);

    gsp->nrs = nRandomStarts;     /* Number of random starts */
    gsp->nts = kNumTempSteps;     /* Number of temperature steps */
    gsp->crs = 50;                /* Number of starts before convergence test */
    gsp->cfm = 20;                /* Top thread frequency convergence criterion */
    gsp->csm = 5;                 /* Top thread start convergence criterion */
    gsp->cet = 5 * SCALING_FACTOR;/* Temperature for convergence test ensemble */
    gsp->cef = 50;                /* Percent of ensemble defining top threads */
    gsp->isl = 1;                 /* Code for choice of starting locations */
    gsp->iso = 0;                 /* Code for choice of segment sample order */
    gsp->ito = 0;                 /* Code for choice of terminus sample order */
    gsp->rsd = -1;                /* Seed for random number generator -- neg for Rand01() */
    gsp->als = 0;                 /* Code for choice of alignment record style */
    gsp->trg = 0;                 /* Code for choice of trajectory record */

    if (fast) {
    //    gsp->nti[0] = 5;            /* Number of iterations per tempeature step */
    //    gsp->nti[1] = 10;
    //    gsp->nti[2] = 25;
        gsp->nti[0] = 10;           /* Number of iterations per tempeature step */
        gsp->nti[1] = 20;
        gsp->nti[2] = 40;
    } else {
        gsp->nti[0] = 10;           /* Number of iterations per tempeature step */
        gsp->nti[1] = 20;
        gsp->nti[2] = 40;
    }

    gsp->nac[0] = 4;              /* Number of alignment cycles per iteration */
    gsp->nac[1] = 4;
    gsp->nac[2] = 4;

    gsp->nlc[0] = 2;              /* Number of location cycles per iteration */
    gsp->nlc[1] = 2;
    gsp->nlc[2] = 2;

    if (fast) {
    //    gsp->tma[0] = 5 * SCALING_FACTOR;  /* Temperature steps for alignment sampling */
    //    gsp->tma[1] = 5 * SCALING_FACTOR;
    //    gsp->tma[2] = 5 * SCALING_FACTOR;
        gsp->tma[0] = 20 * SCALING_FACTOR;  /* Temperature steps for alignment sampling */
        gsp->tma[1] = 10 * SCALING_FACTOR;
        gsp->tma[2] =  5 * SCALING_FACTOR;
    } else {
        gsp->tma[0] = 20 * SCALING_FACTOR;  /* Temperature steps for alignment sampling */
        gsp->tma[1] = 10 * SCALING_FACTOR;
        gsp->tma[2] =  5 * SCALING_FACTOR;
    }

    gsp->tml[0] = 5 * SCALING_FACTOR;     /* Temperature steps for location sampling */
    gsp->tml[1] = 5 * SCALING_FACTOR;
    gsp->tml[2] = 5 * SCALING_FACTOR;

    gsp->lms[0] = 0;              /* Iterations before local minimum test */
    gsp->lms[1] = 10;
    gsp->lms[2] = 20;

    gsp->lmw[0] = 0;              /* Iterations in local min test interval */
    gsp->lmw[1] = 10;
    gsp->lmw[2] = 10;

    gsp->lmf[0] = 0;              /* Percent of top score indicating local min */
    gsp->lmf[1] = 80;
    gsp->lmf[2] = 95;

    if (gsp->als == 0) {
        NumTrajectoryPoints = 1;
    } else {
        NumTrajectoryPoints = 0;
        for (int i=0; i<kNumTempSteps; ++i) {
            NumTrajectoryPoints += gsp->nti[i] * (gsp->nac[i] + gsp->nlc[i]);
        }
        NumTrajectoryPoints *= gsp->nrs;
        NumTrajectoryPoints += gsp->nrs;
    }
    gsp->ntp = NumTrajectoryPoints;

    return(gsp);
}

#define NO_VIRTUAL_COORDINATE(coord) \
    do { coord->type = Threader::MISSING_COORDINATE; return; } while (0)

static void GetVirtualResidue(const AtomSet *atomSet, const Molecule *mol,
    const Residue *res, Threader::VirtualCoordinate *coord)
{
    // find coordinates of key atoms
    const AtomCoord *C = NULL, *CA = NULL, *CB = NULL, *N = NULL;
    Residue::AtomInfoMap::const_iterator a, ae = res->GetAtomInfos().end();
    for (a=res->GetAtomInfos().begin(); a!=ae; ++a) {
        AtomPntr ap(mol->id, res->id, a->first);
        if (a->second->atomicNumber == 6) {
            if (a->second->code == " C  ")
                C = atomSet->GetAtom(ap, true, true);
            else if (a->second->code == " CA ")
                CA = atomSet->GetAtom(ap, true, true);
            else if (a->second->code == " CB ")
                CB = atomSet->GetAtom(ap, true, true);
        } else if (a->second->atomicNumber == 7 && a->second->code == " N  ")
            N = atomSet->GetAtom(ap, true, true);
        if (C && CA && CB && N) break;
    }
    if (!C || !CA || !N) NO_VIRTUAL_COORDINATE(coord);

    // find direction of real or idealized C-beta
    Vector toCB;

    // if C-beta present, vector is in its direction
    if (CB) {
        toCB = CB->site - CA->site;
    }

    // ... else need to calculate a C-beta direction (not C-beta position!)
    else {
        Vector CaN, CaC, cross, bisect;
        CaN = N->site - CA->site;
        CaC = C->site - CA->site;
        // for a true bisector, these vectors should be normalized! but they aren't in other
        // versions of the threader (Cn3D/C and S), so the average is used instead...
//        CaN.normalize();
//        CaC.normalize();
        bisect = CaN + CaC;
        bisect.normalize();
        cross = vector_cross(CaN, CaC);
        cross.normalize();
        toCB = 0.816497 * cross - 0.57735 * bisect;
    }

    // virtual C-beta location is 2.4 A away from C-alpha in the C-beta direction
    toCB.normalize();
    coord->coord = CA->site + 2.4 * toCB;
    coord->type = Threader::VIRTUAL_RESIDUE;

    // is this disulfide-bound?
    Molecule::DisulfideMap::const_iterator ds = mol->disulfideMap.find(res->id);
    coord->disulfideWith =
        (ds == mol->disulfideMap.end()) ? -1 :
        (ds->second - 1) * 2;   // calculate virtualCoordinate index from other residueID
}

static void GetVirtualPeptide(const AtomSet *atomSet, const Molecule *mol,
    const Residue *res1, const Residue *res2, Threader::VirtualCoordinate *coord)
{
    if (res1->alphaID == Residue::NO_ALPHA_ID || res2->alphaID == Residue::NO_ALPHA_ID)
        NO_VIRTUAL_COORDINATE(coord);

    AtomPntr ap1(mol->id, res1->id, res1->alphaID), ap2(mol->id, res2->id, res2->alphaID);
    const AtomCoord
        *atom1 = atomSet->GetAtom(ap1, true, true),   // 'true' means just use first alt coord
        *atom2 = atomSet->GetAtom(ap2, true, true);
    if (!atom1 || !atom2) NO_VIRTUAL_COORDINATE(coord);

    coord->coord = (atom1->site + atom2->site) / 2;
    coord->type = Threader::VIRTUAL_PEPTIDE;
    coord->disulfideWith = -1;
}

static void GetVirtualCoordinates(const Molecule *mol, const AtomSet *atomSet,
    Threader::VirtualCoordinateList *virtualCoordinates)
{
    virtualCoordinates->resize(2 * mol->residues.size() - 1);
    Molecule::ResidueMap::const_iterator r, re = mol->residues.end();
    const Residue *prevResidue = NULL;
    int i = 0;
    for (r=mol->residues.begin(); r!=re; ++r) {
        if (prevResidue)
            GetVirtualPeptide(atomSet, mol,
                prevResidue, r->second, &((*virtualCoordinates)[i++]));
        prevResidue = r->second;
        GetVirtualResidue(atomSet, mol,
            r->second, &((*virtualCoordinates)[i++]));
    }
}

static const int MAX_DISTANCE_BIN = 5;
static int BinDistance(const Vector& p1, const Vector& p2)
{
    double dist = (p2 - p1).length();
    int bin;

    if (dist > 10.0)
        bin = MAX_DISTANCE_BIN + 1;
    else if (dist > 9.0)
        bin = 5;
    else if (dist > 8.0)
        bin = 4;
    else if (dist > 7.0)
        bin = 3;
    else if (dist > 6.0)
        bin = 2;
    else if (dist > 5.0)
        bin = 1;
    else
        bin = 0;

    return bin;
}

static void GetContacts(const Threader::VirtualCoordinateList& coords,
    Threader::ContactList *resResContacts, Threader::ContactList *resPepContacts)
{
    int i, j, bin;

    // loop i through whole chain, just to report all missing coords
    for (i=0; i<coords.size(); ++i) {
        if (coords[i].type == Threader::MISSING_COORDINATE) {
            WARNINGMSG("Threader::CreateFldMtf() - unable to determine virtual coordinate for "
                << ((i%2 == 0) ? "sidechain " : "peptide ") << (i/2));
            continue;
        }

        for (j=i+10; j<coords.size(); ++j) {    // must be at least 10 virtual bonds away

            if (coords[j].type == Threader::MISSING_COORDINATE ||
                // not interested in peptide-peptide contacts
                (coords[i].type == Threader::VIRTUAL_PEPTIDE &&
                 coords[j].type == Threader::VIRTUAL_PEPTIDE) ||
                // don't include disulfide-bonded cysteine pairs
                (coords[i].disulfideWith == j || coords[j].disulfideWith == i)
                ) continue;

            bin = BinDistance(coords[i].coord, coords[j].coord);
            if (bin <= MAX_DISTANCE_BIN) {
                // add residue-residue contact - res1 is lower-numbered residue
                if (coords[i].type == Threader::VIRTUAL_RESIDUE &&
                    coords[j].type == Threader::VIRTUAL_RESIDUE) {
                    resResContacts->resize(resResContacts->size() + 1);
                    resResContacts->back().vc1 = i;
                    resResContacts->back().vc2 = j;
                    resResContacts->back().distanceBin = bin;
                }
                // add residue-peptide contact
                else {
                    resPepContacts->resize(resPepContacts->size() + 1);
                    resPepContacts->back().distanceBin = bin;
                    // peptide must go in vc2
                    if (coords[i].type == Threader::VIRTUAL_RESIDUE) {
                        resPepContacts->back().vc1 = i;
                        resPepContacts->back().vc2 = j;
                    } else {
                        resPepContacts->back().vc2 = i;
                        resPepContacts->back().vc1 = j;
                    }
                }
            }
        }