User talk:Remig

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:

#   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)
}

Contributors

Remig