diff --git a/SCRIPTING.md b/SCRIPTING.md index 2d55526b3..6dac315b2 100644 --- a/SCRIPTING.md +++ b/SCRIPTING.md @@ -181,11 +181,22 @@ BWCENTER ``` If the current working strand is a seven-helix protein (7HP), with Ballesteros-Weinstein numbering for all its transmembrane helices, -then `BWCENTER` will obtain the center of all n.50:CA atom locations and recenter the entire protein. +then `BWCENTER` will obtain the center of all n.50:CA atom locations and recenter the entire strand. This is useful for comparing various active and inactive states of GPCRs; the 1.50 through 7.50 residues will be assumed to be the most stationary parts of the protein, and variations in structure can be observed with minimal global transformational anomalies. +# BWCOPY +Example: +``` +BWCOPY A B +``` + +Copies the transmembrane helix regions and Ballesteros-Weinstein numbers from one strand to another. This is useful if working with two models of +the same protein, one of which is missing the region definitions, for example if one of the models is included with PrimaryDock and the other is not. +The first parameter is the source strand; the second is the destination. + + # CANMOVE Example: ``` diff --git a/data/taar_activations.pepd b/data/taar_activations.pepd new file mode 100644 index 000000000..67965306b --- /dev/null +++ b/data/taar_activations.pepd @@ -0,0 +1,51 @@ + +LOAD "pdbs/TAAR/TAAR1.upright.pdb" A I +UPRIGHT I +BWCENTER + +DOWNLOAD RCSB 8JLN "pdbs/8jln.pdb" ONCE +LOAD "pdbs/8jln.pdb" R A +BWCOPY I A +UPRIGHT A +BWCENTER +SAVE "pdbs/TAAR1.3-iodothyronamine.pdb" + + +DOWNLOAD RCSB 8JLO "pdbs/8jlo.pdb" ONCE +LOAD "pdbs/8jlo.pdb" R A +BWCOPY I A +UPRIGHT A +BWCENTER +SAVE "pdbs/TAAR1.ulotaront.pdb" + + +DOWNLOAD RCSB 8JLP "pdbs/8jlp.pdb" ONCE +LOAD "pdbs/8jlp.pdb" R A +BWCOPY I A +UPRIGHT A +BWCENTER +SAVE "pdbs/TAAR1.ralmitaront.pdb" + + +DOWNLOAD RCSB 8JLQ "pdbs/8jlq.pdb" ONCE +LOAD "pdbs/8jlq.pdb" R A +BWCOPY I A +UPRIGHT A +BWCENTER +SAVE "pdbs/TAAR1.fenoldopam.pdb" + + +DOWNLOAD RCSB 8JLR "pdbs/8jlr.pdb" ONCE +LOAD "pdbs/8jlr.pdb" R A +BWCOPY I A +UPRIGHT A +BWCENTER +SAVE "pdbs/TAAR1.A77636.pdb" + + +DOWNLOAD RCSB 8JSO "pdbs/8jso.pdb" ONCE +LOAD "pdbs/8jso.pdb" R A +BWCOPY I A +UPRIGHT A +BWCENTER +SAVE "pdbs/TAAR1.amphetamine.pdb" diff --git a/src/classes/aminoacid.cpp b/src/classes/aminoacid.cpp index 6bb681f2f..1dc56fa1b 100755 --- a/src/classes/aminoacid.cpp +++ b/src/classes/aminoacid.cpp @@ -2171,23 +2171,25 @@ float AminoAcid::get_omega() return find_angle_along_vector(CA0->get_location(), CA1->get_location(), N->get_location(), axis); } -bond_rotation_fail_reason AminoAcid::rotate_phi(float a) +bond_rotation_fail_reason AminoAcid::rotate_phi(float a, bool absolute) { if (aadef->proline_like) return bf_disallowed_rotation; Atom* N = get_atom("N"); Atom* CA = get_atom("CA"); Bond* b = CA->get_bond_between(N); if (!b) return bf_bond_not_found; + if (absolute) a -= get_phi(); b->rotate(a, true, true); return b->last_fail; } -bond_rotation_fail_reason AminoAcid::rotate_psi(float a) +bond_rotation_fail_reason AminoAcid::rotate_psi(float a, bool absolute) { Atom* CA = get_atom("CA"); Atom* C = get_atom("C"); Bond* b = CA->get_bond_between(C); if (!b) return bf_bond_not_found; + if (absolute) a -= get_psi(); b->rotate(a, true, true); return b->last_fail; } diff --git a/src/classes/aminoacid.h b/src/classes/aminoacid.h index c97575658..aaea71942 100755 --- a/src/classes/aminoacid.h +++ b/src/classes/aminoacid.h @@ -148,8 +148,8 @@ class AminoAcid : public Molecule float get_phi(); float get_psi(); float get_omega(); - bond_rotation_fail_reason rotate_phi(float rel_angle); - bond_rotation_fail_reason rotate_psi(float rel_angle); + bond_rotation_fail_reason rotate_phi(float rel_angle, bool absolute = false); + bond_rotation_fail_reason rotate_psi(float rel_angle, bool absolute = false); // Intermol functions. float get_intermol_binding(AminoAcid* neighbor, bool backbone_atoms_only = false); diff --git a/src/classes/atom.cpp b/src/classes/atom.cpp index 59f6ce583..448304927 100755 --- a/src/classes/atom.cpp +++ b/src/classes/atom.cpp @@ -572,7 +572,7 @@ float Atom::get_charge() { if (Z == 1) { - if (bonded_to && bonded_to[0].btom) + if (bonded_to[0].btom) { float bchg = bonded_to[0].btom->charge; if (!bchg) bchg = bonded_to[0].btom->is_conjugated_to_charge(); diff --git a/src/classes/constants.h b/src/classes/constants.h index a679c7a8b..90fd9fce5 100755 --- a/src/classes/constants.h +++ b/src/classes/constants.h @@ -70,6 +70,7 @@ #define _ALLOW_FLEX_RINGS 0 #define _shield_angle (130.0 * fiftyseventh) #define _shield_angle_pi (100.0 * fiftyseventh) +#define _shield_angle_opposite_charge (65.0 * fiftyseventh) #define _can_clash_angle (180.0 * fiftyseventh) #define _fullrot_stepdeg 20 #define _fullrot_steprad (fiftyseventh*_fullrot_stepdeg) @@ -318,6 +319,11 @@ #define homology_long_axis_rotations 0 #define homology_region_optimization 0 +#define reconnect_length_limit 5 +#define reconnect_generations 10 +#define reconnect_keepbest 5 +#define reconnect_angle_importance 0.1 + ////////////////////////////////////////////////////////////////////////////////////////////// // Debugging stuff. diff --git a/src/classes/molecule.cpp b/src/classes/molecule.cpp index a9622f7c3..7e4c2d052 100644 --- a/src/classes/molecule.cpp +++ b/src/classes/molecule.cpp @@ -423,6 +423,7 @@ void Molecule::hydrogenate(bool steric_only) float valence = atoms[i]->get_valence(); if (valence > 4) valence = 8 - valence; + if (valence > 1) valence += atoms[i]->get_charge(); // cout << atoms[i]->name << " has valence " << valence << endl; @@ -461,7 +462,7 @@ void Molecule::hydrogenate(bool steric_only) char hname[5]; sprintf(hname, "H%d", atcount+1); Atom* H = add_atom("H", hname, atoms[i], 1); - //cout << "Adding " << hname << " to " << atoms[i]->name << " whose valence is " << valence << " and has " << bcardsum << " bonds already." << endl; + // if (!is_residue()) cout << "Adding " << hname << " to " << atoms[i]->name << " whose valence is " << valence << " and has " << bcardsum << " bonds already." << endl; /*atoms[i]->clear_geometry_cache(); SCoord v = atoms[i]->get_next_free_geometry(1); @@ -879,10 +880,16 @@ int Molecule::from_pdb(FILE* is, bool het_only) } esym[1] &= 0x5f; - Point aloc(atof(words[5]), atof(words[6]),atof(words[7])); + int offset = 0; + if (words[4][0] >= 'A' || words[5] < &buffer[29]) offset++; + + Point aloc(atof(words[5+offset]), atof(words[6+offset]),atof(words[7+offset])); + + // cout << esym << " " << aloc << endl; Atom* a = add_atom(esym, words[2], &aloc, 0, 0, charge); added++; + a->pdbidx = atoi(words[1]); // a->residue = atoi(words[4]); @@ -903,6 +910,20 @@ int Molecule::from_pdb(FILE* is, bool het_only) ; } } + else if (!strcmp(words[0], "CONECT")) + { + int i, i1 = atoi(words[1]), i2 = atoi(words[2]); + Atom *a = nullptr, *b = nullptr; + + for (i=0; atoms[i]; i++) + { + if (atoms[i]->pdbidx == i1) a = atoms[i]; + if (atoms[i]->pdbidx == i2) b = atoms[i]; + if (a && b) break; + } + + if (a && b) a->bond_to(b, 1); // TODO: Detect coplanar atoms and add pi bonds where found. + } } buffer[0] = 0; @@ -2192,9 +2213,11 @@ bool Molecule::shielded(Atom* a, Atom* b) const a->shielding_angle = b->shielding_angle = 0; Point aloc = a->get_location(), bloc = b->get_location(); + float achg = a->get_charge(), bchg = b->get_charge(); for (i=0; atoms[i]; i++) { Atom* ai = atoms[i]; + float aichg = ai->get_charge(); if (!ai) break; if (ai == a || ai == b) continue; float rai = ai->distance_to(a); @@ -2207,15 +2230,11 @@ bool Molecule::shielded(Atom* a, Atom* b) const if (f3da > a->shielding_angle) a->shielding_angle = b->shielding_angle = f3da; if (f3da > _shield_angle) { - if (last_iter && (a->residue == 114 || b->residue == 114) && ((a->residue + b->residue) == 114)) - { - /*cout << ai->name << " shields " - << a->residue << ":" << a->name << "..." - << b->residue << ":" << b->name - << " angle " << (f3da*fiftyseven) - << endl;*/ - return true; - } + return true; + } + if (sgn(achg) == sgn(bchg) && sgn(achg) == -sgn(aichg) && f3da > _shield_angle_opposite_charge) + { + return true; } } @@ -2402,7 +2421,7 @@ float Molecule::get_intermol_binding(Molecule** ligands, bool subtract_clashes) ligands[l]->atcount = j; break; } - if (atoms[i]->is_backbone && ligands[l]->atoms[j]->is_backbone + if (atoms[i]->is_backbone && atoms[i]->residue && ligands[l]->atoms[j]->is_backbone && ligands[l]->atoms[j]->residue && ( ( atoms[i]->residue == ligands[l]->atoms[j]->residue - 1 && @@ -2421,14 +2440,20 @@ float Molecule::get_intermol_binding(Molecule** ligands, bool subtract_clashes) float r = ligands[l]->atoms[j]->get_location().get_3d_distance(&aloc); if (r < _INTERA_R_CUTOFF) { + // cout << "🐝 "; if ( !shielded(atoms[i], ligands[l]->atoms[j]) && !ligands[l]->shielded(atoms[i], ligands[l]->atoms[j]) ) { missed_connection.r = 0; + // if (ligands[l]->atoms[j]->residue) cout << ligands[l]->atoms[j]->residue << ":"; + // cout << ligands[l]->atoms[j]->name << " " << ligands[l]->atoms[j]->get_charge() << " ~ "; + // if (atoms[i]->residue) cout << atoms[i]->residue << ":"; + // cout << atoms[i]->name << " " << atoms[i]->get_charge() << " "; // cout << ligands[l]->atoms[j]->get_location().subtract(atoms[i]->get_location()) << ": "; float abind = InteratomicForce::total_binding(atoms[i], ligands[l]->atoms[j]); + // cout << abind << endl; if (abind && !isnan(abind) && !isinf(abind)) { if (abind > 0 && minimum_searching_aniso && ligands[l]->priority) abind *= 1.5; diff --git a/src/classes/point.cpp b/src/classes/point.cpp index ad766d631..b341e0d69 100755 --- a/src/classes/point.cpp +++ b/src/classes/point.cpp @@ -301,7 +301,7 @@ float find_angle_along_vector(Point* pt1, Point* pt2, Point* source, SCoord* v) Point lpt2 = pt2->subtract(source); // Rotate points so v becomes Z axis. - Point cen; + Point cen(0,0,0); Point z(0,0,1); Rotation rots = align_points_3d(&vp, &z, &cen); Point npt1 = rotate3D(&lpt1, &cen, &rots.v, rots.a); diff --git a/src/classes/protein.cpp b/src/classes/protein.cpp index 0e030f167..56e6d9799 100755 --- a/src/classes/protein.cpp +++ b/src/classes/protein.cpp @@ -993,6 +993,13 @@ int Protein::get_bw50(int helixno) return Ballesteros_Weinstein[helixno]; } +void Protein::set_bw50(int hxno, int resno) +{ + if (hxno < 1 || hxno > 78) return; + if (Ballesteros_Weinstein[hxno]) return; + Ballesteros_Weinstein[hxno] = resno; +} + int Protein::get_seq_length() { if (!sequence) return 0; @@ -1554,7 +1561,7 @@ void Protein::rotate_backbone(int resno, bb_rot_dir dir, float angle) } } -void Protein::rotate_backbone_partial(int startres, int endres, bb_rot_dir dir, float angle) +void Protein::rotate_backbone_partial(int startres, int endres, bb_rot_dir dir, float angle, bool bbao) { save_undo_state(); if (startres == endres) return; @@ -1568,8 +1575,50 @@ void Protein::rotate_backbone_partial(int startres, int endres, bb_rot_dir dir, AminoAcid* bendy = get_residue(startres); if (!bendy) return; - LocatedVector lv = bendy->rotate_backbone(dir, angle); - bendy->ensure_pi_atoms_coplanar(); + LocatedVector lv; + if (bbao) + { + Atom *N = bendy->get_atom("N"), *H = bendy->HN_or_substitute(), *CA = bendy->get_atom("CA"), *C = bendy->get_atom("C"), *O = bendy->get_atom("O"); + Bond* b; + switch(dir) + { + case N_asc: + b = N->get_bond_between(CA); + lv = (SCoord)CA->get_location().subtract(N->get_location()); + lv.origin = N->get_location(); + b->rotate(angle, true, true); + break; + + case CA_desc: + b = CA->get_bond_between(N); + lv = (SCoord)CA->get_location().subtract(N->get_location()); + lv.origin = CA->get_location(); + b->rotate(angle, true, true); + break; + + case CA_asc: + b = CA->get_bond_between(C); + lv = (SCoord)CA->get_location().subtract(C->get_location()); + lv.origin = CA->get_location(); + b->rotate(angle, true, true); + break; + + case C_desc: + b = C->get_bond_between(CA); + lv = (SCoord)C->get_location().subtract(CA->get_location()); + lv.origin = C->get_location(); + b->rotate(angle, true, true); + break; + + default: + ; + } + } + else + { + lv = bendy->rotate_backbone(dir, angle); + bendy->ensure_pi_atoms_coplanar(); + } if (lv.r) { @@ -1579,12 +1628,12 @@ void Protein::rotate_backbone_partial(int startres, int endres, bb_rot_dir dir, for (i=startres+inc; movable = get_residue(i); i+=inc) { movable->rotate(lv, angle); - movable->ensure_pi_atoms_coplanar(); + if (!bbao) movable->ensure_pi_atoms_coplanar(); if (i == endres) break; } } - set_clashables(); + if (!bbao) set_clashables(); } void Protein::conform_backbone(int startres, int endres, int iters, bool backbone_atoms_only) @@ -1815,168 +1864,178 @@ void Protein::conform_backbone(int startres, int endres, mass_undoable = wmu; } -#define DBG_BCKCONN 0 -#define _INCREMENTAL_BKCONN 1 -void Protein::backconnect(int startres, int endres) +float Protein::reconnect(int startres, int endres) { - save_undo_state(); - bool wmu = mass_undoable; - mass_undoable = true; - - int i; - int inc = sgn(endres-startres); - - #if _INCREMENTAL_BKCONN - int iter; - for (iter=24; iter>=0; iter--) + int strandlength = abs(endres - startres) + 1, length = strandlength*2; + if (strandlength > reconnect_length_limit) { - #endif - - // Glue the last residue onto the target. - // Then adjust its inner bonds so the other end points as closely to the previous residue as possible. - // Then do the same for the previous residue, and the one before, etc, - // all the way back to the starting residue. - // Give a warning if the starting residue has an anomaly > 0.1A. - AminoAcid *next, *curr, *prev; - int pointer = endres; - float movfactor = 1, decrement, anomaly = 0; - - decrement = 1.0 / fabs(endres - startres); // if exterior is the opposite of interior, what's the opposite of increment? + cout << "Attempt to reconnect strand in excess of length limit." << endl; + throw 0xb5712a9d; + } - next = get_residue(endres+inc); - curr = get_residue(pointer); - prev = get_residue(pointer-inc); + int dir = sgn(endres - startres); if (!dir) dir = 1; - #if DBG_BCKCONN - cout << "backconnect( " << startres << ", " << endres << ")" << endl; - #endif - while (next && curr) - { - #if DBG_BCKCONN - cout << pointer << ":"; - #endif + int connres = (dir<0) ? (endres - 1) : (endres + 1); + AminoAcid* connaa = get_residue(connres); + if (!connaa) + { + cout << "Attempt to reconnect strand to nonexistent target." << endl; + throw 0xb5712a9d; + } - Point pts[4]; - (inc > 0) ? next->predict_previous_COCA(pts) : next->predict_next_NHCA(pts); - Point ptsc[4]; + AminoAcid* endaa = get_residue(endres); + if (!endaa) + { + cout << "Attempt to reconnect discontinuous strand." << endl; + throw 0xb5712a9d; + } - #if _INCREMENTAL_BKCONN - ptsc[0] = (inc > 0) ? curr->get_atom_location("C") : curr->get_atom_location("N"); - ptsc[1] = (inc > 0) ? curr->get_atom_location("O") : curr->HN_or_substitute_location(); - ptsc[2] = curr->get_atom_location("CA"); + char working_candidate[length+1]; + float angle = triangular; + int gen; + int max_candidates = pow(3, length); + float candidates[reconnect_keepbest][max_candidates+1][length+1]; + float best_candidates[reconnect_keepbest][length+1]; + float best_score[reconnect_keepbest]; + float anomaly; - Point _4avg[3]; - for (i=0; i<3; i++) - { - _4avg[0] = pts[i]; - _4avg[1] = ptsc[i]; - _4avg[0].weight = movfactor; - _4avg[1].weight = 1.0 - movfactor; - pts[i] = average_of_points(_4avg, 2); - } - #endif + int ikept, icand, ires, i, j, l; - curr->attach_to_prediction(pts, inc > 0); - // break; + Point should_be_at[4]; + Atom* is_at[4]; + if (dir > 0) + { + connaa->predict_previous_COCA(should_be_at); + is_at[0] = endaa->get_atom("C"); + is_at[1] = endaa->get_atom("O"); + is_at[2] = endaa->get_atom("CA"); + is_at[3] = nullptr; + } + else + { + connaa->predict_next_NHCA(should_be_at); + is_at[0] = endaa->get_atom("N"); + is_at[1] = endaa->HN_or_substitute(); + is_at[2] = endaa->get_atom("CA"); + is_at[3] = nullptr; + } - #if DBG_BCKCONN - cout << "g"; - #endif + for (gen=0; genmovability; - curr->movability = MOV_ALL; - - #if DBG_BCKCONN - cout << "a"; - #endif - - pts[4]; - (inc < 0) ? prev->predict_previous_COCA(pts) : prev->predict_next_NHCA(pts); - Point target_heavy = pts[0]; - Point target_pole = pts[1]; - - #if DBG_BCKCONN - cout << "b"; - #endif + cout << "gen " << gen << " kept " << ikept << " iter " << icand << " "; - float theta, step = fiftyseventh*1.0, r, btheta = 0, bestr; - for (theta=0; theta < M_PI*2; theta += step) + for (ires = 0; ires < length; ires++) { - r = target_heavy.get_3d_distance( (inc > 0) ? curr->get_atom_location("N") : curr->get_atom_location("C") ); - // r -= target_pole.get_3d_distance( (inc > 0) ? curr->HN_or_substitute_location() : curr->get_atom_location("O") ); - if (inc > 0) r -= target_pole.get_3d_distance(curr->HN_or_substitute_location()); - - if (!theta || (r < bestr)) + int resno = startres + dir*(ires>>1); + AminoAcid* aa = get_residue(resno); + if (!aa) { - bestr = r; - btheta = theta; + cout << "Attempt to reconnect discontinuous strand." << endl; + throw 0xb5712a9d; } - curr->rotate_backbone( (inc > 0) ? CA_desc : N_asc, step ); + candidates[ikept][icand][ires] = + ( + gen + ? best_candidates[ikept][ires] // Persist previous best candidates. + : ((ires & 0x1) ? aa->get_psi() : aa->get_phi()) + ) + + angle * working_candidate[ires]; + + cout << candidates[ikept][icand][ires]*fiftyseven << " "; + + // Rotate the backbone. + MovabilityType was_mov = aa->movability; + aa->movability = MOV_ALL; + if (ires & 0x1) rotate_backbone_partial(resno, endres, + (dir>0) ? CA_asc : C_desc, + (candidates[ikept][icand][ires] - aa->get_psi()) * dir, true); + else rotate_backbone_partial(resno, endres, + (dir>0) ? N_asc : CA_desc, + (candidates[ikept][icand][ires] - aa->get_phi()) * dir, true); + aa->movability = was_mov; } - curr->rotate_backbone( (inc > 0) ? CA_desc : N_asc, btheta ); - #if DBG_BCKCONN - cout << "c"; - #endif + // TODO: Score the candidate's anomaly using distance_to_connres + reconnect_angle_importance * anomaly_angle. + anomaly = is_at[0]->get_location().get_3d_distance(should_be_at[0]) + + reconnect_angle_importance * is_at[1]->get_location().get_3d_distance(should_be_at[1]) + + reconnect_angle_importance * is_at[2]->get_location().get_3d_distance(should_be_at[2]) + ; + + cout << anomaly; - btheta=0; - for (theta=0; theta < M_PI*2; theta += step) + // Keep the best reconnect_keepbest from each generation. + for (i=0; i < reconnect_keepbest; i++) { - r = target_heavy.get_3d_distance( (inc > 0) ? curr->get_atom_location("N") : curr->get_atom_location("C") ); - // r -= target_pole.get_3d_distance( (inc > 0) ? curr->HN_or_substitute_location() : curr->get_atom_location("O") ); - if (inc < 0) r -= target_pole.get_3d_distance(curr->get_atom_location("O")); - - if (!theta || (r < bestr)) + if (best_score[i] < 0 || best_score[i] > anomaly) { - bestr = r; - btheta = theta; - } - - curr->rotate_backbone( (inc > 0) ? C_desc : CA_asc, step ); - } - curr->rotate_backbone( (inc > 0) ? C_desc : CA_asc, btheta ); - anomaly = bestr; + for (j=reconnect_keepbest-2; j>=i; j--) + { + best_score[j+1] = best_score[j]; + for (l=0; lmovability = fmov; - } + break; + } + } + + for (ires = 0; ires < length; ires++) + { + working_candidate[ires]++; + if (working_candidate[ires] > 1) working_candidate[ires] = -1; + else break; + } - if (pointer == startres) - { - #if DBG_BCKCONN - cout << ". "; - #endif - break; + cout << endl; } + } - pointer -= inc; - next = curr; - curr = prev; - prev = get_residue(pointer-inc); - movfactor -= decrement; + cout << endl << "Best from generation:" << endl; + for (i=0; i < reconnect_keepbest; i++) + { + cout << i+1 << ": "; + for (l=0; l>1); + AminoAcid* aa = get_residue(resno); - /*if (anomaly > 0.1) cout << "Warning! conform_backbone( " << startres << ", " << endres << " ) anomaly out of range." << endl - << "# " << (startres-inc) << " anomaly: " << anomaly << endl;*/ - #if _INCREMENTAL_BKCONN + // Rotate the backbone. + MovabilityType was_mov = aa->movability; + aa->movability = MOV_ALL; + if (ires & 0x1) rotate_backbone_partial(resno, endres, + (dir>0) ? CA_asc : C_desc, + (best_candidates[0][ires] - aa->get_psi()) * dir, true); + else rotate_backbone_partial(resno, endres, + (dir>0) ? N_asc : CA_desc, + (best_candidates[0][ires] - aa->get_phi()) * dir, true); + aa->movability = was_mov; } - #endif - mass_undoable = wmu; + set_clashables(); + return anomaly; } void Protein::make_helix(int startres, int endres, float phi, float psi) diff --git a/src/classes/protein.h b/src/classes/protein.h index 67bf84de8..5d90b155c 100755 --- a/src/classes/protein.h +++ b/src/classes/protein.h @@ -61,6 +61,7 @@ class Protein void delete_sidechains(int startres, int endres); MetalCoord* coordinate_metal(Atom* metal, int residues, int* resnos, std::vector res_anames); void set_region(std::string name, int start, int end); + void set_bw50(int helixno, int resno); void renumber_residues(int startres, int endres, int new_startres); bool disulfide_bond(int resno1, int resno2); @@ -143,8 +144,8 @@ class Protein LocRotation rotate_piece(int start_res, int end_res, Point origin, SCoord axis, float theta); void rotate_backbone(int residue_no, bb_rot_dir direction, float angle); + void rotate_backbone_partial(int startres, int endres, bb_rot_dir direction, float angle, bool backbone_atoms_only = false); void conform_backbone(int startres, int endres, Atom* a, Point target, int iters = 50); - void rotate_backbone_partial(int startres, int endres, bb_rot_dir direction, float angle); void conform_backbone(int startres, int endres, int iters = 50, bool backbone_atoms_only = false); void conform_backbone(int startres, int endres, Atom* a1, Point target1, Atom* a2, Point target2, Atom* a3, Point target3, int iters = 50); void conform_backbone(int startres, int endres, Atom* a1, Point target1, Atom* a2, Point target2, int iters = 50, bool backbone_atoms_only = false); @@ -155,7 +156,7 @@ class Protein int iters, bool backbone_atoms_only ); - void backconnect(int startres, int endres); + float reconnect(int startres, int endres); void find_residue_initial_bindings(); void undo(); diff --git a/src/classes/scoring.cpp b/src/classes/scoring.cpp index 909b8bd02..549a5e590 100644 --- a/src/classes/scoring.cpp +++ b/src/classes/scoring.cpp @@ -228,7 +228,7 @@ DockResult::DockResult(Protein* protein, Molecule* ligand, Point size, int* addl #endif sprintf(metrics[metcount], "%s%d", reaches_spheroid[i]->get_3letter(), resno); - // cout << metrics[metcount] << ": " << lb << " . "; + // cout << metrics[metcount] << ": " << lb << ". "; if (differential_dock) { diff --git a/src/interpreter.cpp b/src/interpreter.cpp index 1864b3d04..7e8192990 100644 --- a/src/interpreter.cpp +++ b/src/interpreter.cpp @@ -964,7 +964,28 @@ int main(int argc, char** argv) SCoord new_center = old_center; new_center.r *= -1; working->move_piece(1, working->get_end_resno(), new_center); - } // BWCENTER + } // BWCENTER + + else if (!strcmp(words[0], "BWCOPY")) + { + if (!words[1] || !words[2]) raise_error("Insufficient parameters given for BWCOPY."); + + char c = words[1][0]; + if (c < 'A' || c > 'Z') raise_error("Invalid source strand given for BWCOPY."); + Protein* source = strands[c-'A']; + + c = words[2][0]; + if (c < 'A' || c > 'Z') raise_error("Invalid destination strand given for BWCOPY."); + Protein* dest = strands[c-'A']; + + for (l=1; l<=8; l++) + { + std::string rgname = (std::string)"TMR" + std::to_string(l); + int sr = source->get_region_start(rgname), er = source->get_region_end(rgname), bw50 = source->get_bw50(l); + if (sr > 0 && er > 0) dest->set_region(rgname, sr, er); + if (bw50 > 0) dest->set_bw50(l, bw50); + } + } // BWCOPY else if (!strcmp(words[0], "CANMOVE")) { @@ -1073,10 +1094,10 @@ int main(int argc, char** argv) sr = interpret_single_int(words[l++]); if (!words[l]) raise_error("Insufficient parameters given for CONNECT."); ct = interpret_single_int(words[l++]); - if (words[l]) iters = interpret_single_int(words[l++]); + // if (words[l]) iters = interpret_single_int(words[l++]); if (words[l]) raise_error("Too many parameters given for CONNECT."); er = ct - sgn(ct-sr); - +/* Atom *a1, *a2, *a3; AminoAcid *cta = working->get_residue(ct), *era = working->get_residue(er); @@ -1099,7 +1120,9 @@ int main(int argc, char** argv) a3 = era->get_atom("CA"); working->conform_backbone(sr, er, a1, pt3[0], a2, pt3[1], iters); - working->backconnect(sr, er); + working->backconnect(sr, er);*/ + + working->reconnect(sr, er); _no_connect: goto _pc_continue; diff --git a/src/ramachandran.cpp b/src/ramachandran.cpp index afd04932b..916b7657a 100644 --- a/src/ramachandran.cpp +++ b/src/ramachandran.cpp @@ -9,11 +9,23 @@ int main(int argc, char** argv) int i; FILE* fp; bool show_colors = true; + bool skip_helices = false; + bool include_glycine = true, include_not_glycine = true; + bool include_proline = true, include_not_proline = true; + bool include_pre_proline = true, include_not_pre_proline = true; for (i=1; iis_alpha_helix()) continue; + if (aa->get_letter() == 'G') + { + if (!include_glycine || !include_not_proline) continue; + } + else if (aa->get_letter() == 'P') + { + if (!include_proline || !include_not_glycine) continue; + } + else if (!include_not_glycine || !include_not_proline) continue; + + if (!include_pre_proline || !include_not_pre_proline) + { + AminoAcid* next = p.get_residue(i+1); + if (next) + { + if (next->get_letter() == 'P') + { + if (!include_pre_proline) continue; + } + else if (!include_not_pre_proline) continue; + } + } + + float phi = aa->get_phi(), psi = aa->get_psi(); phi = fmod(phi, (M_PI*2)); if (phi >= M_PI) phi -= M_PI*2; diff --git a/src/score_pdb.cpp b/src/score_pdb.cpp index 8451dd31d..a52929756 100644 --- a/src/score_pdb.cpp +++ b/src/score_pdb.cpp @@ -10,6 +10,8 @@ using namespace std; int main(int argc, char** argv) { + _INTERA_R_CUTOFF = _DEFAULT_INTERA_R_CUTOFF; + if (argc < 2) { cout << "Usage:" << endl << endl; @@ -34,13 +36,12 @@ int main(int argc, char** argv) fseek(fp, 0, 0); m.from_pdb(fp, true); - AminoAcid* aa = p.get_residue_bw(3, 50); - if (!aa) aa = p.get_residue_bw(6, 48); - if (!aa) aa = p.get_residue(100); + AminoAcid* aa = p.get_residue(100); if (!aa) aa = p.get_residue(p.get_end_resno()/2); if (aa && aa->get_hydrogen_count() < 3) { + cout << "Hydrogenating..." << endl; int i, n = p.get_end_resno(); for (i=0; i