|
|
| (3 intermediate revisions by the same user not shown) |
| Line 1: |
Line 1: |
| − | I use Jmol to study protein folding. Here is a script I wrote that accepts an amino acid sequence (1 letter encoding: "AAC...FYW" for example) and generates an alpha helix using the Model Kit: | + | What I had here before belongs on the user page - so I moved it... |
| − | <pre># RIBOZOME - Jmol script by Ron Mignery with help from Dr. Angel Herráez
| |
| − | # v1.2 beta 10/23/2013
| |
| − | #
| |
| − | # New in v1.2 beta:
| |
| − | # No longer zaps existing atoms
| |
| − | # Now adds on to existing helix on subsequent runs
| |
| − | # Fixes proline cross-link bug with XXXXPXP
| |
| − | # Localizes variables and globalizes constants
| |
| − | #
| |
| − | # RIBOZOME takes a string message encoding an amino acid (aa) sequence
| |
| − | # and generates a corresponding alpha helix one aa at a time from the
| |
| − | # N terminus to the C terminus rotating the emerging helix as it goes.
| |
| − | #
| |
| − | # The message is a string entered by the user at a prompt.
| |
| − | # It may be typed in or pasted in and be of any length
| |
| − | # If the message is prepended with <C>: (where C is any single letter)
| |
| − | # then the chain is so labeled and separated from existing chains
| |
| − | # if different from the first chain.
| |
| − | #
| |
| − | # The IUPAC/IUBMB 1 letter code is used:
| |
| − | # A=ALAnine B=GLutam?X* C=CYSteine D=ASPartate E=GLUtamate
| |
| − | # F=PHEnylalanine G=GLYcine H=HIStidine I=IsoLEucine K=LYSine
| |
| − | # L=LEUcine M=METhionine N=ASparagiNe O=PYrroLysine*** P=PROline
| |
| − | # Q=GLutamiNe R=ARGinine S=SERine T=THReonine U=SElenoCysteine
| |
| − | # V=VALine W=TRyPtophan X=UNKnown Y=TYRosine Z=ASpar?X**
| |
| − | # *Either GLU or GLN: drawn as GLN with chi3 flipped
| |
| − | # **Either ASP or ASN: drawn as ASN with chi3 flipped
| |
| − | # ***Not supported: drawn as ALA
| |
| − | | |
| − | # The following constant values determine the pitch of the alpha helix
| |
| − | gPHI = -57 # Dihedral angle of N-CA bond (nominally -57)
| |
| − | gPSI = -47 # Dihedral angle of CA-C bond (nominally -47)
| |
| − | gOMEGA = 180 # Dihedral angle of the peptide bond (nominally 180)
| |
| − | gPEPTIDE_ANGLE = 110 # C-N-CA angle (nominally 110)
| |
| − | gPRO_BUMP = -10 # Psi angle change increment to N-3psi when N is Pro
| |
| − | gCHAIN = 'A' # The chain id
| |
| − | | |
| − | # Lookup 3 letter code from 1 letter code
| |
| − | g3from1 = {"A":"ALA", "B":"GLX","C":"CYS", "D":"ASP","E":"GLU", "F":"PHE",
| |
| − | "G":"GLY", "H":"HIS","I":"ILE", "K":"LYS","L":"LEU", "M":"MET",
| |
| − | "N":"GLN", "O":"PYL","P":"PRO", "Q":"GLN","R":"ARG", "S":"SER",
| |
| − | "T":"THR", "U":"SEC","V":"VAL", "W":"TRP","X":"UNK", "Y":"TYR", "Z":"ASX"}
| |
| − | | |
| − | # Generate PDB atom record
| |
| − | function genAtom(e, aa, i, xyz) {
| |
| − | gA = format("ATOM %5d %4s %3s ", gN, e, aa )
| |
| − | gA += format("%s%4d %8.3f", gCHAIN, i, xyz[1] )
| |
| − | gA += format("%8.3f%8.3f\n", xyz[2], xyz[3] )
| |
| − | gN++
| |
| − | return gA
| |
| − | };
| |
| − | | |
| − | # Generate a PDB amino acid record set
| |
| − | function genAA(i, aa) { # Writes globals gA and gN
| |
| − | | |
| − | # From constructed AAs
| |
| − | var N0 = [0.0, 0.0, 0.0]
| |
| − | var CA = [ 0.200, 1.174, 0.911 ]
| |
| − | var C = [ -1.110, 1.668, 1.425 ]
| |
| − | var O = [ -1.320, 1.693, 2.62 ]
| |
| − | var CB = [ 1.062, 2.1950, 0.230 ]
| |
| − |
| |
| − | var G1 = [ 2.359, 1.607, -0.344]
| |
| − | var G2 = [ 1.363, 3.336, 1.157 ]
| |
| − | var D1 = [ 3.222, 2.656, -1.048 ]
| |
| − | var D2 = [ 3.143, 0.904, 0.725 ]
| |
| − | var E1 = [ 3.645, 3.749, -0.167 ]
| |
| − | var E2 = [ 2.491, 3.234, -2.249 ]
| |
| − | var Z = [ 4.470, 4.717, -0.885 ]
| |
| − | var H1 = [ 4.450, 6.006, -0.220 ]
| |
| − | var H2 = [5.833, 4.228, -0.984 ]
| |
| − |
| |
| − | var Gp = [ 2.008, 1.24, -0.46 ]
| |
| − | var Dp = [1.022, 0.213, -1.031 ]
| |
| − |
| |
| − | var Gfy = [ 2.368, 1.471, -0.0152 ]
| |
| − | var D1fy = [ 3.346, 1.524, 0.921 ]
| |
| − | var D2fy = [ 2.493, 0.516, -1.151 ]
| |
| − | var E1fy = [ 4.513, 0.615, 0.8244 ]
| |
| − | var E2fy = [ 3.528, -0.336, -1.206 ]
| |
| − | var Zfy = [ 4.588, -0.285, -0.168 ]
| |
| − | var Hfy = [ 5.738, -1.245, -0.233 ]
| |
| − |
| |
| − | var Ghw = [ 2.406, 1.626, -0.134 ]
| |
| − | var D1hw = [3.498, 1.936, 0.675]
| |
| − | var D2hw = [ 2.713, 0.901, -1.211 ]
| |
| − | var E1hw = [ 4.160, 0.518, -1.178 ]
| |
| − | var E2hw = [ 4.622, 1.160, 0.0816 ]
| |
| − | var E3hw = [ 3.789, 2.523, 1.944 ]
| |
| − | var Z2hw = [ 5.973, 1.177, 0.689 ]
| |
| − | var Z3hw = [ 5.014, 2.550, 2.503 ]
| |
| − | var H2hw = [ 6.153, 1.846, 1.844 ]
| |
| − |
| |
| − | #N1 = [ 2.069, -2.122, -0.554]
| |
| − |
| |
| − | # Build PDB atom records common to all AAs
| |
| − | a3 = g3from1[aa]
| |
| − | if (a3 = "") {
| |
| − | a3 = "UNK"
| |
| − | }
| |
| − | print format("+ %s%d/%d", a3, i, gSeq.count + gResno)
| |
| − | gA = genAtom(" N ", a3, i, N0)
| |
| − | gA += genAtom(" CA ", a3, i, CA)
| |
| − | gA += genAtom(" C ", a3, i, C)
| |
| − | gA += genAtom(" O ", a3, i, O)
| |
| − | if ((aa != 'G') && (aa != 'X')) {
| |
| − | gA += genAtom(" CB ", a3, i, CB)
| |
| − | }
| |
| − | | |
| − | # Now add AA specific atom records
| |
| − | switch (aa) {
| |
| − | case 'A' :
| |
| − | break;
| |
| − | case 'B' :
| |
| − | gA += genAtom(" CG ", a3, i, G1)
| |
| − | gA += genAtom(" CD ", a3, i, D1)
| |
| − | gA += genAtom(" OE1", a3, i, E2) # GLN with Es switched
| |
| − | gA += genAtom(" NE2", a3, i, E1)
| |
| − | break;
| |
| − | case 'C' :
| |
| − | gA += genAtom(" SG ", a3, i, G2)
| |
| − | break;
| |
| − | case 'D' :
| |
| − | gA += genAtom(" CG ", a3, i, G1)
| |
| − | gA += genAtom(" OD1", a3, i, D1)
| |
| − | gA += genAtom(" OD2", a3, i, D2)
| |
| − | break;
| |
| − | case 'E' :
| |
| − | gA += genAtom(" CG ", a3, i, G1)
| |
| − | gA += genAtom(" CD ", a3, i, D1)
| |
| − | gA += genAtom(" OE1", a3, i, E1)
| |
| − | gA += genAtom(" OE2", a3, i, E2)
| |
| − | break;
| |
| − | case 'F' :
| |
| − | gA += genAtom(" CG ", a3, i, Gfy)
| |
| − | gA += genAtom(" CD1", a3, i, D1fy)
| |
| − | gA += genAtom(" CD2", a3, i, D2fy)
| |
| − | gA += genAtom(" CE1", a3, i, E1fy)
| |
| − | gA += genAtom(" CE2", a3, i, E2fy)
| |
| − | gA += genAtom(" CZ ", a3, i, Zfy)
| |
| − | break;
| |
| − | case 'G' :
| |
| − | break;
| |
| − | case 'H' :
| |
| − | gA += genAtom(" CG ", a3, i, Ghw)
| |
| − | gA += genAtom(" ND1", a3, i, D1hw)
| |
| − | gA += genAtom(" CD2", a3, i, D2hw)
| |
| − | gA += genAtom(" CE1", a3, i, E2hw)
| |
| − | gA += genAtom(" NE2", a3, i, E1hw)
| |
| − | break;
| |
| − | case 'I' :
| |
| − | gA += genAtom(" CG1", a3, i, G1)
| |
| − | gA += genAtom(" CG2", a3, i, G2)
| |
| − | gA += genAtom(" CD1", a3, i, D1)
| |
| − | break;
| |
| − | case 'K' :
| |
| − | gA += genAtom(" CG ", a3, i, G1)
| |
| − | gA += genAtom(" CD ", a3, i, D1)
| |
| − | gA += genAtom(" CE ", a3, i, E1)
| |
| − | gA += genAtom(" NZ ", a3, i, Z)
| |
| − | break;
| |
| − | case 'L' :
| |
| − | gA += genAtom(" CG1", a3, i, G1)
| |
| − | gA += genAtom(" CD1", a3, i, D1)
| |
| − | gA += genAtom(" CD2", a3, i, D2)
| |
| − | break;
| |
| − | case 'M' :
| |
| − | gA += genAtom(" CG ", a3, i, G1)
| |
| − | gA += genAtom(" SD ", a3, i, D1)
| |
| − | gA += genAtom(" CE ", a3, i, E1)
| |
| − | break;
| |
| − | case 'N' :
| |
| − | gA += genAtom(" CG ", a3, i, G1)
| |
| − | gA += genAtom(" OD1", a3, i, D1)
| |
| − | gA += genAtom(" ND2", a3, i, D2)
| |
| − | break;
| |
| − | case 'P' :
| |
| − | gA += genAtom(" CG ", a3, i, GP)
| |
| − | gA += genAtom(" CD ", a3, i, DP)
| |
| − | break;
| |
| − | case 'Q' :
| |
| − | gA += genAtom(" CG ", a3, i, G1)
| |
| − | gA += genAtom(" CD ", a3, i, D1)
| |
| − | gA += genAtom(" OE1", a3, i, E1)
| |
| − | gA += genAtom(" NE2", a3, i, E2)
| |
| − | break;
| |
| − | case 'R' :
| |
| − | gA += genAtom(" CG ", a3, i, G1)
| |
| − | gA += genAtom(" CD ", a3, i, D1)
| |
| − | gA += genAtom(" NE ", a3, i, E1)
| |
| − | gA += genAtom(" CZ ", a3, i, Z)
| |
| − | gA += genAtom(" NH1", a3, i, H1)
| |
| − | gA += genAtom(" NH2", a3, i, H2)
| |
| − | break;
| |
| − | case 'S' :
| |
| − | gA += genAtom(" OG ", a3, i, G1)
| |
| − | break;
| |
| − | case 'T' :
| |
| − | gA += genAtom(" OG1", a3, i, G1)
| |
| − | gA += genAtom(" CG2", a3, i, G2)
| |
| − | break;
| |
| − | case 'U' :
| |
| − | gA += genAtom("SeG ", a3, i, G1)
| |
| − | break;
| |
| − | case 'V' :
| |
| − | gA += genAtom(" CG1", a3, i, G1)
| |
| − | gA += genAtom(" CG2", a3, i, G2)
| |
| − | break;
| |
| − | case 'W' :
| |
| − | gA += genAtom(" CG ", a3, i, Ghw)
| |
| − | gA += genAtom(" CD1", a3, i, D1hw)
| |
| − | gA += genAtom(" CD2", a3, i, D2hw)
| |
| − | gA += genAtom(" CE2", a3, i, E2hw)
| |
| − | gA += genAtom(" NE1", a3, i, E1hw)
| |
| − | gA += genAtom(" CE3", a3, i, E3hw)
| |
| − | gA += genAtom(" CZ2", a3, i, Z2hw)
| |
| − | gA += genAtom(" CZ3", a3, i, Z3hw)
| |
| − | gA += genAtom(" CH2", a3, i, H2hw)
| |
| − | break;
| |
| − | case 'X' :
| |
| − | gA += genAtom(" Xx ", a3, i, CB)
| |
| − | break;
| |
| − | case 'Y' :
| |
| − | gA += genAtom(" CG ", a3, i, Gfy)
| |
| − | gA += genAtom(" CD1", a3, i, D1fy)
| |
| − | gA += genAtom(" CD2", a3, i, D2fy)
| |
| − | gA += genAtom(" CE1", a3, i, E1fy)
| |
| − | gA += genAtom(" CE2", a3, i, E2fy)
| |
| − | gA += genAtom(" CZ ", a3, i, Zfy)
| |
| − | gA += genAtom(" OH ", a3, i, Hfy)
| |
| − | break;
| |
| − | case 'Z' :
| |
| − | gA += genAtom(" CG ", a3, i, G1)
| |
| − | gA += genAtom(" OD1", a3, i, D2) # ASN with Ds switched
| |
| − | gA += genAtom(" ND2", a3, i, D1)
| |
| − | break;
| |
| − | default :
| |
| − | break;
| |
| − | }
| |
| − | | |
| − | return gA
| |
| − | };
| |
| − | | |
| − | # Generate an alpha helix
| |
| − | function genAlpha(gSeq) {
| |
| − | | |
| − | gN = all.count + 1 # global new N atom index
| |
| − | | |
| − | # Find last linkable N if any
| |
| − | gResno = 0 # global pre-existing AA count
| |
| − | var pn = 1 # previous gN
| |
| − | for (var i = all.count-1; i > 0; i--) {
| |
| − |
| |
| − | # If found
| |
| − | if (distance({atomno=i}, {0,0,0}) < 0.1) {
| |
| − | pn = i
| |
| − |
| |
| − | # If new chain, separate from existing chain
| |
| − | if ({atomno=i}.chain != gCHAIN) {
| |
| − | select all
| |
| − | translateselected {2.069, -2.122, -0.554 } #N1
| |
| − | }
| |
| − | else {
| |
| − | gResno = {atomno=i}.resno
| |
| − | }
| |
| − | break;
| |
| − | }
| |
| − | }
| |
| − | | |
| − | # For each aa
| |
| − | set appendnew false
| |
| − | var nn = gN # new N
| |
| − | for (var i = 1; i <= gSeq.count; i++) {
| |
| − | | |
| − | # Move polypeptide C to bond distance from new AA N
| |
| − | select all
| |
| − | fix none
| |
| − | translateselected {2.069, -2.122, -0.554 } #N1
| |
| − | | |
| − | # Gen AA
| |
| − | gA = "data \"append aa\"\n" # global PDB atom record
| |
| − | gA += genAA(i + gResno, gSeq[i]); # gN is updated in subroutine
| |
| − | gA += "end \"append aa\""
| |
| − | script inline @{gA}
| |
| − | | |
| − | # If PRO
| |
| − | var pb = 0
| |
| − | if ((gSeq.count - i) >= 2) {
| |
| − | if (gSeq[i + 2] == 'P') {
| |
| − | pb = gPRO_BUMP
| |
| − | }
| |
| − | }
| |
| − | | |
| − | # If not first AA
| |
| − | if (nn > 1) {
| |
| − | | |
| − | # Gen axis on new N perpendicular to the plane
| |
| − | # containing atoms nn, pn+2 and nn+1
| |
| − | var v1={atomno = @{pn+2}}.xyz - {atomno = nn}.xyz
| |
| − | var v2={atomno = @{nn+1}}.xyz - {atomno = nn}.xyz
| |
| − | var axis = cross(v1, v2)
| |
| − |
| |
| − | # Center on atom previous C
| |
| − | axis += {atomno = @{pn+2}}.xyz
| |
| − | | |
| − | # Rotate the polypeptide N on the new AA C to tetrahedral (nominally 110)
| |
| − | select atomno < nn
| |
| − | fix atomno >= nn
| |
| − | rotateselected @axis {atomno = nn} @{gPEPTIDE_ANGLE - 69.3}
| |
| − |
| |
| − | # Make omega dihedral = gOMEGA (nominally 180)
| |
| − | rotateselected {atomno=@{pn+2}} {atomno=nn} @{gOMEGA - 147.4}
| |
| − | | |
| − | # Make the new phi dihedral = gPHI (nominally -57)
| |
| − | rotateselected {atomno = nn} {atomno = @{nn+1}} @{gPHI - 8.7}
| |
| − | | |
| − | # Make the old psi dihedral = gPSI (nominally -47)
| |
| − | select atomno < @{pn+2} && atomno != @{pn+3}
| |
| − | rotateselected {atomno=@{pn+1}} {atomno=@{pn+2}} @{gPSI + 33.4 + pb}
| |
| − | }
| |
| − |
| |
| − | # Step new and previous N
| |
| − | pn = nn
| |
| − | nn = gN
| |
| − | | |
| − | # Make the peptide bond
| |
| − | connect
| |
| − | }
| |
| − |
| |
| − | # Clean up
| |
| − | connect ([UNK].CA) ([UNK].Xx and within(group, _1))
| |
| − | select all
| |
| − | fix none
| |
| − | print format("%d atoms generated", gN)
| |
| − | }
| |
| − | | |
| − | echo Generating Alpha Helix
| |
| − | | |
| − | # Get the sequence from the user
| |
| − | gSeq = prompt("Enter AA sequence (1 char coded)", "")%9999%0
| |
| − | if (gSeq.count > 0) {
| |
| − | if (gSeq[2] == ':') {
| |
| − | gCHAIN = gSeq[1]
| |
| − | gSeq[1] = ' '
| |
| − | gSeq[2] = ' '
| |
| − | gSeq = gSeq%0
| |
| − | }
| |
| − | print format ("Sequence=%s phi=%d psi=%d", gSeq, gPHI, gPSI)
| |
| − | print format ("chain=%s peptide angle=%d pro bump=%d", gCHAIN, gPEPTIDE_ANGLE, gPRO_BUMP)
| |
| − | genAlpha(gSeq)
| |
| − | }
| |
| − | </pre>
| |
