Difference between revisions of "User talk:Remig"

From Jmol
Jump to navigation Jump to search
m
(Replaced content with "What I had here before belongs on the user page - so I moved it...")
 
(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>
 

Latest revision as of 21:32, 31 October 2013

What I had here before belongs on the user page - so I moved it...

Contributors

Remig