The Genetic Code and DNA Geometry

Alex Kasman
College of Charleston

Abstract

This Mathematica notebook is designed to accompany the article of the same name.  By downloading the notebook, a user can either look at and test my code to verify the claims made in the paper or can use the commands for their own research.

Note that the computation of the Mutual Information in Section 5 includes a very time consuming step.  To avoid reproducing that step, one can separately download the values that it generates in the file “numberinboxgivencodon.txt.zip”, unzip and load it instead to simply set the values directly.

Section 1: Introduction

1.1 Duplexed Codes (Some Examples)

Using the function “exampf” I can send either a 0 or a 1 to a recipient by sending them a number from 1 to 4 with the understanding they will apply this command to it.

exampf[s_]:=Which[s==1,0,s==2,0,s==3,1,s==4,1]

Similarly, I can send either a -1 or a 1 to the a recipient by sending them a 1, 2, 3, or 4 if I know they will apply this command “exampg” to it:

exampg[s_]:=Which[s==1,-1,s==2,1,s==3,-1,s==4,1]

The real point of this example is that I can use the same “coded message” s to send any message I want to both recipients.  For example, if I want to send a 1 to the first recipient and a -1 to the second, I can send them both “s=3” and that will be decoded in the way I want.

exampf[3]

exampg[3]

On the other hand, these codes are not duplexed as well:

exampF[s_]:=Which[s==1,1,s==2,0,s==3,1,s==4,1]
exampG[s_]:=Which[s==1,-1,s==2,1,s==3,-1,s==4,1]

Again, I can send either one of them a message that I want, but I can no longer send both of them an arbitrary message using the same s.  In particular, there is no input for which the output of the first function is 0 and the output of the second function is -1.

We will later see that “mutual information” is one way to quantify the effectiveness of duplexing, but the purpose here is just to introduce the concept.

1.2 The Genetic Code

There are four possible bases in a DNA sequence.  A sequence of length 2 is called a dimer and a sequence of length 3 is called a codon.

bases={"A","C","G","T"}

dimers=Tuples[bases,2]

codons=Tuples[bases,3]

Now, I define two functions: one that gives the index number of a given codon in that sequence and the other, an inverse function of it, that gives the sequence associated to an index.

codon[i_]:=codons[[i]]
codonnumof[codon_]:=Position[codons,codon][[1,1]]

codon[3]

codonnumof[{"A","A","G"}]

Breaking a codon down into its constituent dimers:

delta[{a_,b_,c_}]:={a,b}
delta'[{a_,b_,c_}]:={b,c}

Let’s test it on:

codon[2]

delta[codon[2]]

delta'[codon[2]]

The genetic code assigns an amino acid (or the word “Stop”, which I will call an amino acid also just for simplicity even though it really is not one) to each codon.  These commands define that correspondence:

setcodon[liszt_,b_]:=Module[{i},
For[i=1,i≤Length[liszt],i=i+1,PhiI[liszt[[i]]]=b];
PhiIInv[b]=liszt;
aminoacids=Append[aminoacids,b];]

aminoacid[i_]:=aminoacids[[i]]
aminoacidnumfor[v_]:=Position[aminoacids,v][[1,1]]

aminoacids={};

setcodon[{{"A","T","T"},{"A","T","C"},{"A","T","A"}},"I"]

setcodon[{{"C","T","T"},{"C","T","C"},{"C","T","A"},{"C","T","G"},{"T","T","A"},{"T","T","G"}},"L"]

setcodon[{{"G","T","T"},{"G","T","C"},{"G","T","A"},{"G","T","G"}},"V"]

setcodon[{{"T","T","T"},{"T","T","C"}},"F"]

setcodon[{{"A","T","G"}},"M"]

setcodon[{{"T","G","T"},{"T","G","C"}},"C"]

setcodon[{{"G","C","T"},{"G","C","C"},{"G","C","A"},{"G","C","G"}},"A"]


setcodon[{{"G","G","T"},{"G","G","C"},{"G","G","A"},{"G","G","G"}},"G"]


setcodon[{{"C","C","T"},{"C","C","C"},{"C","C","A"},{"C","C","G"}},"P"]


setcodon[{{"A","C","T"},{"A","C","C"},{"A","C","A"},{"A","C","G"}},"T"]


setcodon[{{"T","C","T"},{"T","C","C"},{"T","C","A"},{"T","C","G"},{"A","G","T"},{"A","G","C"}},"S"]


setcodon[{{"T","A","T"},{"T","A","C"}},"Y"]


setcodon[{{"T","G","G"}},"W"]


setcodon[{{"C","A","A"},{"C","A","G"}},"Q"]


setcodon[{{"A","A","T"},{"A","A","C"}},"N"]


setcodon[{{"C","A","T"},{"C","A","C"}},"H"]


setcodon[{{"G","A","A"},{"G","A","G"}},"E"]


setcodon[{{"G","A","T"},{"G","A","C"}},"D"]


setcodon[{{"A","A","A"},{"A","A","G"}},"K"]

setcodon[{{"C","G","T"},{"C","G","C"},{"C","G","A"},{"C","G","G"},{"A","G","A"},{"A","G","G"}},"R"]


setcodon[{{"T","A","A"},{"T","A","G"},{"T","G","A"}},"Stop"]

The commands above have defined several things for us.  First, there is the set of amino acids:

aminoacids

Then, a function from the set of codons to that called PhiI which is exactly the famous “genetic code”

PhiI[codon[3]]

That is true because K is the amino acid corresponding to the codon AAG.  It is common to view this correspondence by listing the pre-image for each of the amino acids:

MatrixForm[Table[{aminoacid[i],PhiIInv[aminoacid[i]]},{i,1,Length[aminoacids]}]]

Section 2: DNA Geometry

2.1 Position Vectors for Rungs

A “rung” in a DNA molecule will be represented by a rectangle in 3-dimensional space whose short and long sides have lengths given below.  In addition, I will draw in a little bit of a vector orthogonal to the rectangle with length veclength:

shortsidelength=10
longsidelength=20
veclength=3

respectively.  (The size of the rectangle doesn’t really matter.  It just needs to be drawn large enough that you can see its orientation and small enough so that they do not intersect much.)

A position vector for a rung will be (x1,x2,x3,x4,x5,x6) where  (x1,x2,x3) is just the usual spacial coordinates (x, y and z) of the center of the rectangle and (x4,x5,x6) are Tait-Bryan angles giving its orientation.   It is useful to “project” onto those separately:

pi1[{x1_,x2_,x3_,x4_,x5_,x6_}]:={x1,x2,x3}
pi2[{x1_,x2_,x3_,x4_,x5_,x6_}]:={x4,x5,x6}

Tait-Bryan Angles

The command TBtoFrame[{x4,x5,x6}] returns {v1,v2,v3} where the vector v1 will be parallel to the short sides, v2 parallel to the long sides and v3 will be orthogonal to it (following the right - hand - rule) and together (v1, v2, v3) form an orthonormal frame.

TBtoFrame[{x4_,x5_,x6_}]:=Transpose[{{Cos[x5] Cos[x6],-Cos[x5] Sin[x6],Sin[x5]},{Cos[x6] Sin[x4] Sin[x5]+Cos[x4] Sin[x6],Cos[x4] Cos[x6]-Sin[x4] Sin[x5] Sin[x6],-Cos[x5] Sin[x4]},{-Cos[x4] Cos[x6] Sin[x5]+Sin[x4] Sin[x6],Cos[x6] Sin[x4]+Cos[x4] Sin[x5] Sin[x6],Cos[x4] Cos[x5]}}]

For small x4, x5 and x6 with x6 non-negative, this can be inverted using:

FrametoTB[{v1_,v2_,v3_}]:=Module[{x4,x5,x6},
x5=ArcSin[v3[[1]]];
x4=ArcSin[-v3[[2]]/Cos[x5]];
x6=ArcCos[v1[[1]]/Cos[x5]];
{x4,x5,x6}]

TBtoFrame[{.1,.2,1.3}]

FrametoTB[%]

Drawing a rung in space:

When I draw a rung in space, I ask which base it is so that I can color different bases differently, but that was not mentioned in the paper.

basecolor["A"]:=RGBColor[1,0,0]
basecolor["T"]:=RGBColor[0,1,0]
basecolor["C"]:=RGBColor[0,0,1]
basecolor["G"]:=RGBColor[1,1,0]

showrungframe[{center_,frame_},base_]:=Module[{v1,v2,v3},
v1=frame[[1]];
v2=frame[[2]];
v3=frame[[3]];
Show[Graphics3D[{basecolor[base],Polygon[{center+shortsidelength/2 v1+ longsidelength/2 v2,center-shortsidelength/2 v1 +longsidelength/2 v2,center-shortsidelength/2 v1 - longsidelength/2 v2, center +shortsidelength/2 v1-longsidelength/2 v2}],AbsoluteThickness[2],GrayLevel[0],Line[{center,center+veclength v3}]}]]]

showrung[p_,base_]:=showrungframe[{pi1[p],TBtoFrame[pi2[p]]},base]

So, for instance, this draws two rungs.  One in the “standard starting position” of {0,0,0,0,0,0} and another that is 5 angstroms up and then tilted by pi/8 in around the x-axis.

Show[showrung[{0,0,0,0,0,0},"A"],showrung[{0,0,5,Pi/8,0,0},"T"]]

2.2 Hassan-Calladine Parameters

This section defines a command DimerStep(p,s) where p=(x1,...,x6) is a position vector for a rung and s is a 6-vector of Hassan-Calladine step parameters.  The output of DimerStep(p,s) is (x1’,...,x6’), the position vector for another rung, the unique rung whose relationship to the first is given by those Hassan-Calladine parameters.  In other words, this command performs the transformation that turns a rung into the next rung as given by the step parameters.

First, I just need the rotation matrices whose arguments are given in degrees:

Rx[theta_]:={{1,0,0},{0,DCos[theta],-DSin[theta]},{0,DSin[theta],DCos[theta]}}

Ry[theta_]:={{DCos[theta],0,DSin[theta]},{0,1,0},{-DSin[theta],0,DCos[theta]}}

Rz[theta_]:={{DCos[theta],-DSin[theta],0},{DSin[theta],DCos[theta],0},{0,0,1}}

DCos[x_]:=N[Cos[Pi x/180]]
DSin[x_]:=N[Sin[Pi x/180]]

The command DimerStep takes a position vector with Tait-Bryan angles as its first argument, but it calls a more fundamental command that is defined using a frame for the orientation.  That way it is possible to avoid the conversion to and from angles which only work for a restricted set of angles.

DimerStepFrame[{oldcenter_,oldframe_},s_]:=Module[{Tmat,Tmidmat,newcenter,Gamma,phi,oldu1,oldu2,midv1,midv2,midnorm,coordchange,coordchangeinv,shift,slide,rise,tilt,roll,twist},
shift=s[[1]];slide=s[[2]];rise=s[[3]];tilt=s[[4]];roll=s[[5]];twist=s[[6]];
coordchange=Transpose[oldframe];
Gamma=Sign[roll]Sqrt[roll^2+tilt^2];
phi=If[Gamma==0,phi=0,phi=ArcSin[tilt/Gamma]*180/Pi];
Tmat=coordchange.Rz[twist/2-phi].Ry[Gamma].Rz[twist/2+phi];
Tmidmat=coordchange.Rz[twist/2-phi].Ry[Gamma/2].Rz[phi];
midv1=Tmidmat.{1,0,0};
midv2=Tmidmat.{0,1,0};
midnorm=Cross[midv1,midv2];
newcenter=oldcenter+shift midv1 + slide midv2 + rise midnorm;
{newcenter,Transpose[Tmat]}]

DimerStep[p_,s_]:=Module[{oldcenter,oldframe,newcenter,newframe,rem},
oldframe=TBtoFrame[pi2[p]];
oldcenter=pi1[p];
rem=DimerStepFrame[{oldcenter,oldframe},s];
newcenter=rem[[1]];
newframe=rem[[2]];
Flatten[{newcenter,FrametoTB[newframe]}]]

So, for example, the rung whose relationship to the rung {0,0,0,0,0,0} has Hassan-Calladine parameters {1,2,3,4,5,6} (i.e. shift=1, slide=2, rise=3, tilt=4, roll=5 and twist=6) is:

DimerStep[{0,0,0,0,0,0},{1,2,3,4,5,6}]

They look like this:

Show[showrung[{0,0,0,0,0,0},"A"],showrung[DimerStep[{0,0,0,0,0,0},{1,2,3,4,5,6}],"T"]]

2.3 Dimer Geometry (the data from Olson et al)

sbar[{"C","G"}]={0,.41,3.39,0,5.4,36.1}
sbar[{"C","A"}]={.09,.53,3.33,.5,4.7,37.3}
sbar[{"T","A"}]={0,.05,3.42,0,3.3,37.8}
sbar[{"A","G"}]={.09,-.25,3.34,-1.7,4.5,31.9}
sbar[{"G","G"}]={.05,-.22,3.42,-.1,3.6,32.9}
sbar[{"A","A"}]={-.03,-.08,3.27,-1.4,.07,35.1}
sbar[{"G","A"}]={-.28,.09,3.37,-1.5,1.9,36.3}
sbar[{"A","T"}]={0,-.59,3.31,0,1.1,29.3}
sbar[{"A","C"}]={.13,-.58,3.36,-.1,.7,31.5}
sbar[{"G","C"}]={0,-.38,3.4,0,.3,33.6}

shat[{"C","G"}]={.87,.56,.27,4.2,5.2,5.5}
shat[{"C","A"}]={.55,.89,.26,3.7,5.1,6.5}
shat[{"T","A"}]={.52,.71,.24,2.7,6.6,5.5}
shat[{"A","G"}]={.69,.41,.23,3.3,3.4,4.5}
shat[{"G","G"}]={.76,.64,.24,3.7,4.5,5.2}
shat[{"A","A"}]={.57,.45,.22,3.3,5.4,3.9}
shat[{"G","A"}]={.46,.70,.26,3.8,5.3,4.4}
shat[{"A","T"}]={.57,.31,.21,2.5,4.9,4.5}
shat[{"A","C"}]={.59,.41,.23,3.1,3.9,4.2}
shat[{"G","C"}]={.61,.56,.24,3.9,4.6,4.7}

The parameters for the six dimers not given in the definitions above can be determined from their "complements":

switchtiltandshift[v_]:={-v[[1]],v[[2]],v[[3]],-v[[4]],v[[5]],v[[6]]}

sbar[{"T","T"}]=switchtiltandshift[sbar[{"A","A"}]]
sbar[{"T","G"}]=switchtiltandshift[sbar[{"A","C"}]]
sbar[{"C","T"}]=switchtiltandshift[sbar[{"G","A"}]]
sbar[{"C","C"}]=switchtiltandshift[sbar[{"G","G"}]]
sbar[{"T","C"}]=switchtiltandshift[sbar[{"A","G"}]]
sbar[{"G","T"}]=switchtiltandshift[sbar[{"C","A"}]]

the standard deviations do not change under this complementarity

shat[{"T","T"}]=shat[{"A","A"}]
shat[{"T","G"}]=shat[{"A","C"}]
shat[{"C","T"}]=shat[{"G","A"}]
shat[{"C","C"}]=shat[{"G","G"}]
shat[{"T","C"}]=shat[{"A","G"}]
shat[{"G","T"}]=shat[{"C","A"}]

Now, this command will perform a dimer step as before, but instead of giving the step parameters, you give a dimer “d” and a vector of “z-scores” that say by how many standard deviations the step parameter exceeds its expected value for that dimer:

gamma[p_,dimer_,zvec_]:=DimerStep[p,sbar[dimer]+DiagonalMatrix[zvec].shat[dimer]]

Here you can see how changing the z-score changes things.  In this sequence, the expected geometry is the one in the middle figure and the others are variations in which some of the z-scores are taking values between -4 and 4.

Table[Show[showrung[{1,2,3,.1,.2,.1},"A"],showrung[gamma[{1,2,3,.1,.2,.1},{"G","C"},{zscore,0,-zscore,zscore,0,0}],"T"]],{zscore,-4,4,1}]

2.4 Expected Sequence Geometry (including Expected Codon Geometry)

In the paper I define “Gamma-Bar” to be a function that gives the expected geometry for a codon (a sequence of length 3) assuming its initial rung is in position {0,0,0,0,0,0}.  Here, I’ll define it to be more generally something that gives as output a 3-vector whose first component is the sequence, second component is the expected position of the last rung (as in the paper) and whose third component is a picture of the molecule in its expected geometry

[Notice that the command only converts to Tait-Bryan angles at the end, using frames until then.  That means that the image should be correct even if the angles become large, but the position vector given at the end may not be correct if the molecule has bent out of range.  That is not a problem for codons but would be for the long bent molecule in the example.]

GammaBar[sequence_]:=Module[{i,currentp,currentbase,dimer,graph},
currentp={{0,0,0},IdentityMatrix[3]};
currentbase=sequence[[1]];
graph=showrungframe[currentp,currentbase];
For[i=2,i≤Length[sequence],i=i+1,
dimer={sequence[[i-1]],sequence[[i]]};
currentbase=sequence[[i]];
currentp=DimerStepFrame[currentp,sbar[dimer]];
graph=Show[graph,showrungframe[currentp,currentbase]];];
{sequence,Flatten[{currentp[[1]],FrametoTB[currentp[[2]]]}],Show[graph,PlotRange->All]}]

A good test sequence for it is a repeat of a bunch of these “bending sequences” because they are known to bend tightly.

bendingsequence[n_]:=Flatten[Table[{"A","A","A","A","A","C","G","G","G","C"},{i,1,n}]]

GammaBar[bendingsequence[7]]

For codons, the expected geometries are:

expectedcodongeomtable=Table[GammaBar[codon[i]],{i,1,64}];

MatrixForm[expectedcodongeomtable]

2.5 Deformable Codon Geometry

The command in 2.4 only gives the EXPECTED geometry based on the average found in the Olson data.   In other words, it completely ignores the fact that the molecule can bend, and that the amount it can bend is also sequence dependent.  This information is contained in the standard deviations from Section 2.3.  So, here we define a command GammaDef that takes not only a codon but also two sets of z-scores as arguments and it returns the geometry for the third rung assuming the Hassan-Calladine parameters exceed their expected values by those z-scores multiplied by the standard deviations:

GammaDef[codon_,zvec_,zprimevec_]:=gamma[gamma[{0,0,0,0,0,0},delta[codon],zvec],delta'[codon],zprimevec]

Notice that the first six z-scores are applied to the first dimer and the next six are applied to the second dimer.  

When the z-scores are all zero, this just gives the expected geometry of the third rung, just like GammaBar would if it was applied to the codon sequence:

GammaDef[{"A","C","T"},{0,0,0,0,0,0},{0,0,0,0,0,0}]

GammaBar[{"A","C","T"}]

But, the z-scores can vary corresponding to bending of the DNA molecule within the ranges found reasonable by Olson et al as reflected in their standard deviations:

GammaDef[{"A","C","T"},{1,2,-.2,3,1,-1},{-2,1,1,-1,.5,.5}]

Section 3: Alternative Genetic Codes / Permutations

Okay, so what I want to do is define a function sigma[] that takes in a codon and spits out another codon, and it should have the property that the first two bases of codon1 and codon2 match iff the first two bases of sigma[codon1] and sigma[codon2] match.

So, basically, there is a permuation first of the 16 PAIRS of bases

It could be useful to have commands that turn the index number i of one of the 64 codons into two numbers: a CLASS that identifies the first two bases and a subclass that identifies the third.

codonsubclassof[codon_]:=Mod[codonnumof[codon]-1,4]
codonclassof[codon_]:=(codonnumof[codon]-codonsubclassof[codon]-1)/4

codonclassof[{"G","C","T"}]

codonsubclassof[{"G","C","T"}]

Table[{codon[i],codonclassof[codon[i]],codonsubclassof[codon[i]]},{i,1,64}]

Because of the order that the codons are listed in, we have a simple inversion 4 codonclassof[codon[i]]+ codonsubclass[codon[i]] + 1 = i:

classtonum[{i_,j_}]:=4 i + j + 1

Table[classtonum[{codonclassof[codon[i]],codonsubclassof[codon[i]]}],{i,1,64}]

Right.

So, I guess that what I want is a permutation first of the 16 classes, and then I want a separate random permutation of the four bases within each class.

What I’ll do is define classpermutation to be a permutation of {0,..,.15} and subclasspermutation[c] to be a permutation of {0,1,2,3} for each c in {0,...,15}.  

The default would be the identity:

permutationisidentity:=(classpermutation=Table[i,{i,0,15}];For[i=0,i≤15,i=i+1,subclasspermutation[i]=Table[j,{j,0,3}]])

permutationisidentity

classpermutation

subclasspermutation[1]

But, if the command pickrandompermutation is called, it resets those to something else:

giveapermutationof[n_]:=Module[{unpicked,i,j,outlist,rand},
unpicked=Table[i,{i,0,n-1}];
outlist={};
For[i=0,i≤n-1,i=i+1,rand=Random[Integer,{1,Length[unpicked]}];
outlist=Append[outlist,unpicked[[rand]]];
unpicked=Delete[unpicked,rand]];
outlist]

giveapermutationof[16]

giveapermutationof[4]

pickrandompermutation:=Module[{i},
classpermutation=giveapermutationof[16];
For[i=0,i≤15,i=i+1,
subclasspermutation[i]=giveapermutationof[4];
]]

Okay, so what I want to do is define a function sigma[] that takes in a codon and spits out another codon, and it should have the property that the first two bases of codon1 and codon2 match iff the first two bases of sigma[codon1] and sigma[codon2] match.

So, basically, there is a permuation first of the 16 PAIRS of bases

It could be useful to have commands that turn the index number i of one of the 64 codons into two numbers: a CLASS that identifies the first two bases and a subclass that identifies the third.

codonsubclassof[codon_]:=Mod[codonnumof[codon]-1,4]
codonclassof[codon_]:=(codonnumof[codon]-codonsubclassof[codon]-1)/4

codonclassof[{"G","C","T"}]

codonsubclassof[{"G","C","T"}]

Table[{codon[i],codonclassof[codon[i]],codonsubclassof[codon[i]]},{i,1,64}]

It should be true that 4 codonclassof[codon[i]]+ codonsubclass[codon[i]] + 1 = i, right

classtonum[{i_,j_}]:=4 i + j + 1

Table[classtonum[{codonclassof[codon[i]],codonsubclassof[codon[i]]}],{i,1,64}]

Right.

So, I guess that what I want is a permutation first of the 16 classes, and then I want a separate random permutation of the four bases within each class.

What I’ll do is define classpermutation to be a permutation of {0,..,.15} and subclasspermutation[c] to be a permutation of {0,1,2,3} for each c in {0,...,15}.  

The default would be the identity:

permutationisidentity:=(classpermutation=Table[i,{i,0,15}];For[i=0,i≤15,i=i+1,subclasspermutation[i]=Table[j,{j,0,3}]])

permutationisidentity

classpermutation

subclasspermutation[1]

But, if the command pickrandompermutation is called, it resets those to something else:

giveapermutationof[n_]:=Module[{unpicked,i,j,outlist,rand},
unpicked=Table[i,{i,0,n-1}];
outlist={};
For[i=0,i≤n-1,i=i+1,rand=Random[Integer,{1,Length[unpicked]}];
outlist=Append[outlist,unpicked[[rand]]];
unpicked=Delete[unpicked,rand]];
outlist]

giveapermutationof[16]

giveapermutationof[4]

pickrandompermutation:=Module[{i},
classpermutation=giveapermutationof[16];
For[i=0,i≤15,i=i+1;
subclasspermutation[i]=giveapermutationof[4];
]]

The permutation “sigma” that turns a codon into a corresponding codon is then:

permuteclass[incodon_]:=classpermutation[[codonclassof[incodon]+1]]

permutesubclass[incodon_]:=subclasspermutation[permuteclass[incodon]][[codonsubclassof[incodon]+1]]

sigma[incodon_]:=codon[classtonum[{permuteclass[incodon],permutesubclass[incodon]}]]

sigma[{"A","C","T"}]

sigma[{"A","C","A"}]

pickrandompermutation

sigma[{"A","C","T"}]

sigma[{"A","C","A"}]

The “alternative genetic code” would then just be sigma followed by the real genetic code:

Phisigma[codon_]:=PhiI[sigma[codon]]

Phisigma[{"A","C","T"}]

Phisigma[{"A","C","A"}]

And, we need an inverse function.  Probably the easiest way to do it is just this:

PhisigmaInv[amino_]:=Module[{rem,i},rem={};
For[i=1,i≤64,i=i+1,If[Phisigma[codon[i]]==amino,rem=Append[rem,codon[i]]]];
rem]

Here is an illustration of its use:

pickrandompermutation

PhisigmaInv["P"]

Phisigma[{"T","A","A"}]

Phisigma[{"T","A","C"}]

pickrandompermutation

PhisigmaInv["P"]

Phisigma[{"A","A","C"}]

Phisigma[{"A","A","T"}]

Section 4: Total Network Length

The graphs we will be considering have these vertices, either the centers of the expected rungs:

V1=Table[pi1[expectedcodongeomtable[[i,2]]],{i,1,Length[codons]}]

or the Tait-Bryan angles of the expected rungs:

V2=Table[pi2[expectedcodongeomtable[[i,2]]],{i,1,Length[codons]}]

maxV1={Max[Table[V1[[i,1]],{i,1,Length[V1]}]],Max[Table[V1[[i,2]],{i,1,Length[V1]}]],Max[Table[V1[[i,3]],{i,1,Length[V1]}]]}

minV1={Min[Table[V1[[i,1]],{i,1,Length[V1]}]],Min[Table[V1[[i,2]],{i,1,Length[V1]}]],Min[Table[V1[[i,3]],{i,1,Length[V1]}]]}

Notice that the third component varies by less than .3 while the other two vary by more than 1.  That is why I project out the third component in the graphs shown in the paper.

The distance formula will just be ordinary Euclidean distance

dist[{v1_,v2_}]:=Sqrt[(v1-v2).(v1-v2)]

This adds up the distances for all pairs in a set:

totaldist[liszt_]:=Module[{allpairs,i},
allpairs=Subsets[liszt,{2}];
Sum[dist[allpairs[[i]]],{i,1,Length[allpairs]}]]

The sets that I will be adding up over will be the vertices that correspond to codons with the same amino acid under the map Phisigma

PhisigmaInv["P"]

totaldistanceforamino[amino_]:=Module[{codonlist,pointlist,i},
codonlist=PhisigmaInv[amino];
pointlist=Table[V1[[codonnumof[codonlist[[i]]]]],{i,1,Length[codonlist]}];
totaldist[pointlist]]

permutationisidentity

totaldistanceforamino[aminoacid[4]]

totaldistanceforallaminos:=Sum[totaldistanceforamino[aminoacid[i]],{i,1,Length[aminoacids]}]

permutationisidentity

totaldistanceforallaminos

rem=Table[pickrandompermutation;totaldistanceforallaminos,{i,1,4000}];

rem

<<HypothesisTesting`

MeanCI[rem]

Mean[rem]

Max[rem]

Min[rem]

totaldistanceforamino2[amino_]:=Module[{codonlist,pointlist,i},
codonlist=PhisigmaInv[amino];
pointlist=Table[V2[[codonnumof[codonlist[[i]]]]],{i,1,Length[codonlist]}];
totaldist[pointlist]]

permutationisidentity

totaldistanceforallaminos2:=Sum[totaldistanceforamino2[aminoacid[i]],{i,1,Length[aminoacids]}]

permutationisidentity

totaldistanceforallaminos2

rem2=Table[pickrandompermutation;totaldistanceforallaminos2,{i,1,4000}];

rem2

<<HypothesisTesting`

MeanCI[rem2]

Mean[rem2]

Max[rem2]

Min[rem2]

4.3 Comparing with a Specific Alternative

The goal here is to generate an example which has more geometric freedom for a given amino acid.  This is what the GPH would have predicted is true for the real genetic code.  Apparently, that’s not the case, and the goal now is just to find an example which demonstrates that it would have been possible to do much “better” than the code we’ve actually got. I’m not looking for THE extreme example, but just AN extreme example .  Hopefully one for which it will be visibly obvious how and why it has more geometric freedom than the real genetic code.

liszt will be a list of pairs that are the totaldistanceforallaminos and totaldistandforallaminos2 for each random permutation

while I’m going along, I’ll keep track of the largest WEIGHTED SUM of the two so far

scalar=Max[rem]/Max[rem2]

Then, whenever a bigger sum comes along, I’ll record that permutation with this command:

storecurrentpermutation[td_,td2_,i_]:=Module[{j},(Print["New max achieved ->",{td,td2},"=",td+scalar td2];
biggestsum=td+scalar td2;
biggestct=i;
extremepermutation={classpermutation,Table[subclasspermutation[j],{j,0,15}]})]

doone[i_]:=Module[{t1,t2},pickrandompermutation;
t1=totaldistanceforallaminos;
t2=totaldistanceforallaminos2;
If[t1+scalar t2>biggestsum,storecurrentpermutation[t1,t2,i]];
liszt=Append[liszt,{t1,t2}];]

For[i=0;biggestct=0;biggestsum=0;liszt={},i≤400,i=i+1,doone[i+1]]

extremepermutation

permutationisextreme:=(classpermutation=extremepermutation[[1]];For[i=0,i≤15,i=i+1,subclasspermutation[i]=extremepermutation[[2,i+1]]];)

permutationisextreme

Now the permutation is set to be that extreme one.  We can check by looking at the total distances and they should be as shown in the last printed line above:

totaldistanceforallaminos

totaldistanceforallaminos2

seetotaldistforaminoacid[i_]:=Module[{codonlist,pointlist,j},
codonlist=PhisigmaInv[aminoacid[i]];
pointlist=Table[V1[[codonnumof[codonlist[[j]]]]],{j,1,Length[codonlist]}];
showpointlist2d[pointlist,RGBColor[1,0,0],{.00,.00,.00}]]

showpointlist2d[rawpointlist_,color_,dashing_]:=Module[{subsets,j,pointlist},
pointlist=Table[{rawpointlist[[j,1]],rawpointlist[[j,2]]},{j,1,Length[rawpointlist]}];
If[Length[pointlist]==1,subsets={{pointlist[[1]],pointlist[[1]]}},
subsets=Subsets[pointlist,{2}]];
Show[Graphics[{AbsoluteThickness[1],color,Table[Line[subsets[[j]]],{j,1,Length[subsets]}]}],ListPlot[pointlist,PlotStyle->AbsolutePointSize[6]],PlotRange->{{-.5,1.0},{-1.1,1.0}},Axes->True]]

permutationisidentity

This is the REAL genetic code.

Show[Table[seetotaldistforaminoacid[i],{i,1,Length[aminoacids]}]]

It has two isolated points that are both near the bottom right:

Show[seetotaldistforaminoacid[5],seetotaldistforaminoacid[13]]

But the extreme permutation looks like this:

permutationisextreme

Show[Table[seetotaldistforaminoacid[i],{i,1,Length[aminoacids]}]]

It still has two isolated points (it MUST because of the permutations we have chosen) but they are now closer to the middle:

Show[seetotaldistforaminoacid[5],seetotaldistforaminoacid[13]]

We can also look at the graphs whose vertices are the Tait-Bryan angles:

points=Show[Graphics3D[Table[{Red,PointSize[.02],Point[V2[[i]]]},{i,1,Length[V2]}]]]

However, this time it is not clear which variable I would project out, and so instead I will actually view it in 3D.

showpointlist3d[pointlist_,color_,dashing_]:=Module[{subsets,j},
subsets=Subsets[pointlist,{2}];
Show[Graphics3D[{AbsoluteThickness[1],color,Table[Line[subsets[[j]]],{j,1,Length[subsets]}]}],PlotRange->All,Axes->False]]

seetotaldistforaminoacid2[i_]:=Module[{codonlist,pointlist,j},
codonlist=PhisigmaInv[aminoacid[i]];
pointlist=Table[V2[[codonnumof[codonlist[[j]]]]],{j,1,Length[codonlist]}];
showpointlist3d[pointlist,RGBColor[1,0,0],{.00,.020,.010}]]

permutationisidentity

Show[points,Table[seetotaldistforaminoacid2[i],{i,1,Length[aminoacids]}]]

permutationisextreme

Show[points,Table[seetotaldistforaminoacid2[i],{i,1,Length[aminoacids]}]]

Again, not only does it just LOOK like there is more length in the lines, you can see that the isolated point at the bottom in the real code is linked to other codon choices in the alternative genetic code.  This all suggests that the real genetic code was NOT selected to be one with particularly wide geometric choice for any given DNA sequence.

Section 5: Mutual Information

Finding the Ranges of the Six Geometric Parameters Over 3 SD Deformations

sdtable=Tuples[{-3,3},6];

maxrem=GammaDef[codon[1],sdtable[[1]],sdtable[[1]]]

minrem=%;

checkmaxmin[newin_]:=Module[{i,tester},tester=False;
For[i=1,i≤6,i=i+1,
If[newin[[i]]>maxrem[[i]],maxrem[[i]]=newin[[i]];tester=True;];
If[newin[[i]]<minrem[[i]],minrem[[i]]=newin[[i]];tester=True];
]
If[tester,Print[MatrixForm[{minrem,maxrem}]]];
]

This will run through all of the codons and all combinations of standard deviations with absolute value 3, keeping track of the largest and smallest value for each of the six geometric parameters:

Do[checkmaxmin[GammaDef[codon[i1],sdtable[[i2]],sdtable[[i3]]]],{i1,1,Length[codons]},{i2,1,Length[sdtable]},{i3,1,Length[sdtable]}]

Then, the parameters go from

minrem

to

maxrem

Putting Positions into Bins

we need 6-dimensional boxes, three for the x, y, z coordinates and three for the angles.  The min and max values are set by the vectors minrem and maxrem computed above.

how many boxes?  (Larger is harder and more accurate...)

nn=4;

now, we want to be able to tell which "box" a given set of coordinates is in

beta[vect_]:=Table[betai[vect[[i]],i],{i,1,6}]
betai[u_,i_]:=Module[{max,min,skip,rem},
max=maxrem[[i]];
min=minrem[[i]];
skip=(max-min)/(nn);
rem=1+Floor[(u-min)/skip];
Min[{nn,Max[{rem,1}]}]] (* the Min and Max at the end are just to make sure it never gives a number less than one or bigger than nn *)

Then, given the position for the third rung of a codon, this converts it into a six-tuple of elements of {1,2,3,4} indicating where this lies relative to the ranges for each.  For instance:

beta[GammaDef[codon[5],{-4,2,1,2,1,4},{-2,3,1,4,2,4}]]

This tells us that the angle x6 for this is near the top of the range, but that x1, x2 and x4 are below the middle of their ranges.

Numerically Generating a Probability Distribution

Ideally, I’d like a function P[omega | C] that gives the probability the geometry ends up in the box with indices omega given that the codon is C.  But, I’m not sure there is an easy analytic way to do that.  Instead, I’m going to generate it numerically by running through lots of possible configurations (varying the two step parameters) and just determining which box it ends up in.  I’ll then step a counter for that box by the probability of that configuration (assuming normal distributions for the Olson parameters).

First, what is that probability factor that comes from the normal approximation?

One Dimension: If we have a segment of width dx that is centered at z standard deviations from the mean the probability will be

NormProb[z_,dx_]:=NIntegrate[1/Sqrt[2Pi]Exp[-x^2/2],{x,z-dx/2,z+dx/2}]

NormProb[0,2]

That’s right..empirical rule says that is the area within one standard deviation of the mean

NormProb[0,4]

2*NormProb[-1,2]

That seems right!

Now, let’s assume we will be looking zn boxes starting with one centered at -zmax and ending at zmax.

zmax=3
zn=4

dz=2zmax/(zn-1)

So, we’re talking about having boxes centered at

zcenters=Table[-zmax+i dz,{i,0,zn-1}]

The ranges they represent are

Table[{-zmax-dz/2+i dz,-zmax +dz/2 +i dz},{i,0,zn-1}]

So, the probabilities we need are:

Clear[prob]
For[i=0,i<zn,i=i+1,prob[-zmax + i dz]=NormProb[-zmax + i dz,dz]]

The percentage of the area under the normal curve we are theoretically addressing then is:

Sum[prob[-zmax + i dz],{i,0,zn-1}]

And, for vectors (representing the z scores in each of the 6 parameters for two dimers) it will just be a product:

prob[vec1_,vec2_]:=Module[{i},Product[prob[vec1[[i]]],{i,1,Length[vec1]}]Product[prob[vec2[[i]]],{i,1,Length[vec2]}]]

Now, I need a count for each box index (integers from 1 to nn for each of six indices) and each codon number (from 1 to 64).

by default, they are all zero:

So, I’m going to loop through the geometries for all nn values for each parameter, note the box, and add the associated probability to that box for that codon:

addboxcount[vec1_,vec2_,cdn_]:=Module[{boxnums,pr},
boxnums=beta[GammaDef[codon[cdn],vec1,vec2]];
pr=prob[vec1,vec2];
numberinboxgivencodon[boxnums,cdn]=numberinboxgivencodon[boxnums,cdn]+pr;
]

The standard deviations I will consider are all of the ones that are in {-3,-1,3,1}:

sdtable=Tuples[zcenters,6];

Since this is too much to compute all at once, I will save the output in files:

DIR="/RESEARCH/DNAGeometry"

For me, the easiest way to do this was to split it between two computers and have Mathematica split it between four processors.  These commands are commented out because they take a very long time to run.  You should download the file numberinboxgivencodon.txt.zip, uncompress it, and execute it to set the values.  But, if you would like to check that it works you can uncomment these and execute them.

(* 
ParallelDo[Clear[numberinboxgivencodon];
numberinboxgivencodon[boxnums_,cdn_]:=0;For[ct2=1,ct2≤Length[sdtable],ct2=ct2+1,For[ct1=1,ct1≤Length[sdtable],ct1=ct1+1,addboxcount[sdtable[[ct1]],sdtable[[ct2]],ct3]]];
Save[DIR<>"/geometryforcodon"<>ToString[ct3]<>".txt",numberinboxgivencodon];
Print["Completed codon ",ct3],{ct3,1,32}] *)

(* 
ParallelDo[Clear[numberinboxgivencodon];
numberinboxgivencodon[boxnums_,cdn_]:=0;For[ct2=1,ct2≤Length[sdtable],ct2=ct2+1,For[ct1=1,ct1≤Length[sdtable],ct1=ct1+1,addboxcount[sdtable[[ct1]],sdtable[[ct2]],ct3]]];
Save[DIR<>"/geometryforcodon"<>ToString[ct3]<>".txt",numberinboxgivencodon];
Print["Up to codon ",ct3],{ct3,33,64}]*)

<<numberinboxgivencodon.txt

Computing the Mutual Information

Some definitions needed for MI

The formulas in Information Theory seem to assume that x Log y is zero when y is zero, but Mathematica does not think so (it considers it undefined).  Hence, I’ll define a command that will handle this properly.

LogProd[x_,y_]:=Which[x==0,x,True,x Log[2,y]]

Then, the formulas for Mutual Information involve sums over the two sets of output values of the codes.  In this case, those will be the set of amino acids

aminoacids

and the set of all “bin indices” which is called Omega in the paper:

Omega=Tuples[{1,2,3,4},6]

Let me define commands that will sum over these.  The way it works is that you define a function FUNC[x] and then say SumOveraminoacids[FUNC] or SumOverOmega[FUNC].

SumOveraminoacids[FUNC_]:=Module[{i},Sum[FUNC[aminoacids[[i]]],{i,1,Length[aminoacids]}]]

SumOverOmega[FUNC_]:=Module[{i},Sum[FUNC[Omega[[i]]],{i,1,Length[Omega]}]]

SumOveraminoacidsandOmega[FUNC_]:=Module[{i,j},Sum[FUNC[aminoacids[[i]],Omega[[j]]],{i,1,Length[aminoacids]},{j,1,Length[Omega]}]]

Here’s an illustration.  The command “Mean” (built into Mathematica) turns an element of Omega into a number.

Mean[{1,2,3,1,4,1}]

Then this command will compute the sum of all of those.

SumOverOmega[Mean]

The Probabilities We Need

The probability of getting a certain amino acid:

Pa[amino_]:=Length[PhisigmaInv[amino]]/64

Although this technically depends on the permutation sigma, it actually doesn’t because the permutations do not change the number of codons for a given amino acid:

permutationisidentity

Pa["M"]

Pa["Y"]

pickrandompermutation

Pa["M"]

Pa["Y"]

So, it is actually better to define it using PhiIInv because that is computationally faster:

Clear[Pa]
Pa[amino_]:=Length[PhiIInv[amino]]/64

No matter what the codon was, the total amount spread out between all of the boxes for that given codon was:

(*totalinallboxesforgivencodon=Sum[prob[sdtable[[i]],sdtable[[j]]],{i,1,Length[sdtable]},{j,1,Length[sdtable]}]*)

totalinallboxesforgivencodon=0.9992401549500121`

It makes sense that this is almost 1 because I was counting all of the probability between -4 and 4 standard deviations so little is outside of it.  Still, I will divide by that when computing the probability to get a properly normalized value for the probability of a certain box given a codon.

Pomegagivenc[omega_,c_]:=numberinboxgivencodon[omega,c]/totalinallboxesforgivencodon

So, as both a test and an illustration, this is how I would add up P(omega|c) over all of omega for a given codon c:

test[omega_]:=Pomegagivenc[omega,5]

SumOverOmega[test]

As expected, that’s 1.

Then P(omega) would be

Pomega[omega_]:=Sum[Pomegagivenc[omega,c],{c,1,64}]/64

Finally, the joint probability function is

Paandomega[a_,omega_]:=Module[{codonlist,i},
codonlist=PhisigmaInv[a];
Sum[Pa[a]Pomegagivenc[omega,codonnumof[codonlist[[i]]]]/Length[codonlist],{i,1,Length[codonlist]}]]

NB: This is the one thing in the computation of MI that actually does depend on the selection of the permutation

permutationisidentity

Paandomega["Y",{2,2,2,2,2,2}]

pickrandompermutation

Paandomega["Y",{2,2,2,2,2,2}]

pickrandompermutation

Paandomega["Y",{2,2,2,2,2,2}]

Computing the Entropies of the Two Codes

For the genetic code the entropy is:

Hasummand[x_]:=LogProd[Pa[x],Pa[x]]

Ha=-N[SumOveraminoacids[Hasummand]]

And for the geometric code it is:

Homegasummand[y_]:=LogProd[Pomega[y],Pomega[y]]

Homega=-N[SumOverOmega[Homegasummand]]

Then, the mutual information of the two codes must be between 0 and

Min[Ha,Homega]

The Mutual Information of the two codes

MIsummand[x_,y_]:=LogProd[Paandomega[x,y],Paandomega[x,y]]-LogProd[Paandomega[x,y],(Pa[x]Pomega[y])]

The mutual information in the case of the REAL genetic code is:

permutationisidentity

realMI=SumOveraminoacidsandOmega[MIsummand]

And, we can do the same for n=50 randomly selected alternative codes to see how this compares:

SetSharedVariable[Liszt]

Liszt={}

ParallelDo[pickrandompermutation; AppendTo[Liszt,SumOverXandY[MIsummand]],50]

Liszt

Min[Liszt]

Mean[Liszt]

Max[Liszt]

MeanCI[Liszt]